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ABSTRACT 


An  ablation  type  constrictor  has  been  developed  in  which  an  electric  arc  at  very  high 
pressure  is  sufficiently  stable  for  spectrographic  measurements.  Radial  temperature  dis¬ 
tributions  were  determined  for  250  ampere  arcs  in  Oelrin,  (COH2)n,  plasma  at  100  and  150 
atmospheres.  Voltage  gradient  is  constant  along  the  arc  axis  and  as  high  as  500  volts/cm. 
Pressure  gradient  in  the  constrictor  is  negligible.  These  two  results  indicate  strongly  that 
in  spite  of  the'  self-induced  flow  field,  the  ablation  type  arc  is  homogeneous  along  the  column. 

T empe'rature  measurement  was  based  on  non-optically  thin  Abel  inversion  of  continuum 
intensity  mea'surements  useful  to  an  optical  depth  of  one  or  more,  which  related  temperature 
to  theoretically  computed  continuum  emission.  A  previously  developed  arc  radiation  model 
was  extended  to  include  a  few  selected  lines,  and  these  contribute  significantly  to  total  radia¬ 
tion  because  they  are  very  much  broadened. 
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SECTION  I 


INTRODUCTION 


One  of  the  critical  problems  to  be  dealt  with  in  the  development  of  high  enthalpy,  high 
pressure  electric  arc  heaters,  and  also  in  the  generation  of  high  pressure  arc  plasmas  for 
basic  studies  of  transport  properties,  is  the  strong  influence  of  radiation  in  the  transfer  of 
energy  from  the  arc  to  the  surrounding  gas  and  to  the  container  walls.  The  objective  of  the 
study  was  to  determine  temperature  profiles  and  heat  fluxes  from  such  high  pressure  arcs 
(100  atmospheres  and  up).  In  the  pressure  range  of  interest,  self-absorption  is  an  important 
phenomenon  not  only  for  discrete  spectral  lines  but  increasingly  even  for  continuum  radiation. 
It  must  therefore  be  accounted  for  in  an  analysis  of  the  heat  flux  from  such  an  arc. 


Under  a  previous  contract,  an  analytical  model  and  computer  program  (called  ARCRAD  III) 
were  developed  for  computation  of  radiant  interchange  within,  and  net  radiant  emission  from, 
a  cylindrically  symmetric  constricted  arc.  The  importance  of  self -absorption  was  shown  to  be 
very  much  dependent  on  frequency.  In  the  vacuum  ultraviolet  part  of  the  spectrum,  no  energy 
reaches  the  constrictor  wall  while  only  a  minor  correction  was  necessary  in  computing  tem¬ 
perature  from  continuum  intensity  measurements  in  the  infrared  (at  8330  fi). 


Radiant  emission  from  high  pressure  arcs  is  diminished  somewhat  by  self-absorption  but 
still  constitutes  a  serious  problem  for  arc  confinement.  Earlier  workW  showed,  for  instance, 
specific  radial  heat  fluxes  up  to  100  lcw/cm2  at  the  constrictor  wall.  This  value  far  exceeds 
the  upper  energy  limit  (~10  kw/cm2)  for  steady  state,  wall-stabilized  arcs  with  water  cooled 
copper  constrictors  so  a  practical  very  high  pressure  arc  apparatus  must  use  either  stabilization 
principles.  The  ablation  type  arc  has  been  used  for  the  highest  energies  and  pressures^2-7). 
Other  successful  experiments  have  been  made  with  a  "free  burning  arc,"(8>9)  with  arcs  in 
rapid  rotation  in  magnetic  fields^10)  and  with  "swirl"  stabilized  arcs. 


Whereas  it  seems  unlikely  that  a  steady  state,  wall-stabilized,  arcing  apparatus  can 
successfully  be  operated  at  very  high  pressures  and  currents,  the  advantages  of  the  wall- 
stabilized  column  for  a  quantitative  measurement  of  the  plasma  parameters  are  obvious. 

Our  early  efforts  ^  were  therefore  directed  toward  transient  operation  of  wall-stabilized 
arcs  in  small  diameter,  uncooled  copper  constrictor  cascades.  Serious  instabilities  of  the 
plasma  column  were  encountered  at  pressures  of  the  order  of  100  atm. ,  particularly  at  high 
currents,  (400  amp),  which  limited  our  measurements  to  a  study  of  the  simpler  arc  para¬ 
meters  such  as  current,  voltage  gradient,  arc  diameter  and  radial  heat  flux. 

A  measurement  of  the  more  important  plasma  parameter  temperature  was  accomplished 
at  a  current  of  100  amps  at  a  strongly  reduced  specific  radiative  flux.  At  this  condition,  the 
arc  was  still  relatively  stable. 


To  expand  on  the  arc  parameters  pressure  and  current,  in  the  work  covered  by  this  report, 
two  approaches  were  followed: 
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First,  working  with  a  copper  cage  constrictor,  additional  effort  was  made  to  eliminate 
disturbances  and  thus  delay  the  onset  of  instabilities;  also,  external  magnetic  fields  were 
applied  to  supplement  the  stabilization  afforded  by  the  constrictor.  Second,  an  ablation  type 
constrictor  was  built.  This  is  a  constrictor  made  of  a  material  which  is  an  electrical  insula¬ 
tor  and  which  evaporates  under  action  of  the  arc.  Of  these,  only  the  ablation  type  constrictor 
proved  to  be  sufficiently  stable  for  measurement  of  the  temperature  distribution. 

Two  wall  materials  were  used  for  the  ablation  type  constrictor:  Plexiglas  for  preliminary 
work  because  its  transparency  permitted  photographic  observation  of  the  entire  arc  column, 
and  Delrin  because  its  composition  is  almost  pure  (H2CO)n  (it  does  not  carbonize),  and  thermo¬ 
dynamic  and  radiative  properties,  while  neither  simple  nor  immediately  available,  offered 
some  prospect  of  yielding  to  calculation. 

Progressing  to  higher  arc  pressure  and  current  and  therefore  to  higher  density  and  higher 
temperature  also  meant  that  the  effect  of  self-absorption  on  intensity  profile  measurements 
would  be  inescapable.  It  therefore  becomes  essential  to  account  for  it  in  the  Abel  inversion 
process  of  converting  intensity  to  temperature  profiles.  A  procedure  suggested  by  Griem^11^ 
was  developed  and  proved  highly  successful. 

Use  of  the  continuum  for  temperature  measurements  required  investigation  of  the 
contribution  to  total  absorption  coefficient  from  the  wings  of  nearby  spectral  lines.  A  few 
lines  were  also  chosen  to  check  their  contribution  to  net  radiant  emission. 
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SECTION  II 


ANALYSIS 

A.  RADIANT  EMISSION  FROM  SE  LF- ABSORBING  ARC  COLUMN 

The  model  used  in  computing  radiant  emission  from  a  self-absorbing  arc  column,  and 
the  analysis  based  on  it,  have  been  described  in  detail^.  The  radiating  arc  column  is 
assumed  to  be  cylindrically  symmetric,  with  a  length-to-diameter  ratio  sufficiently  large 
(10  or  more)  to  neglect  end  effects.  The  arc  cylinder  is  divided  into  annuli,  each  at  a 
uniform  temperature  and  density.  Then,  taking  into  account  the  effect  of  self-absorption 
within  the  annulus,  volumetric  emission  from  each  annulus  is  converted  to  an  equivalent 
radiant  flux  from  the  annulus  boundary  surfaces.  The  assumption  of  local  thermodynamic 
equilibrium  (LTE)  permits  the  application  of  Kirchoff's  law  relating  absorption,  transmission, 
and  emission.  Except  for  the  core,  each  annulus  has  an  interior  (concave)  surface  as  well  as 
an  exterior  (convex)  surface  from  which  radiant  emission  must  be  taken  into  account. 

Radiant  interchange  among  all  annuli  is  computed  by  considering  emission,  absorption, 
and  transmission  along  a  representative  array  of  paths  from  all  the  emitting  surfaces  to  the 
outside  of  the  column.  Because  self-absorption  is  significant,  the  controlling  parameter  is 
optical  depth  which  is  the  dimensionless  product  of  absorption  coefficient  and  a  characteristic 
length.  Two  optical  depths  are  distinguished  here:  1)  absorption  lengths  along  the  trans¬ 
mission  path,  and  2)  emission  lengths  within  each  annulus. 

Computation  of  radiant  interchange  among  the  annuli  proceeded  as  follows,  Figure  1. 
Consider  a  surface  element  dA  of  an  annulus,  radiating  with  an  average  intensity  I  within  a 
solid  angle  do>  over  a  wave  number  interval  n  in  a  direction  making  an  angle  y  with  the  normal 
to  dA.  This  equivalent  surface  emission  arises  from  the  volume  within  the  annuli  contained 
within  the  projection  of  do;  back  through  dA.  Some  of  the  radiated  energy  will  be  absorbed  by 
other  annuli  along  the  path  and  some  will  escape  the  arc  column  entirely.  The  radiant  inter¬ 
change  between  every  possible  pair  of  annuli  was  computed,  taking  into  account  the  optical 
depth  of  the  absorption  path  between  the  annuli.  After  appropriate  summations,  which  are  a 
massive  sorting  and  summing  job  ideally  suited  for  a  digital  computer,  the  result  is  net 
emission  for  each  annulus.  The  fact  that  not  all  transmission  paths  intersect  all  annuli 
results  in  the  absence  of  some  terms,  but  these  were  taken  care  of  in  the  scheme  for  indexing 
the  summation  loops  of  computer  program  ARCRAD  mb. 

This  program  is  essentially  the  same  as  ARCRAD  III,  which  was  presented  in  detail  in 
the  previous  report,  except  that: 

(1)  To  meet  the  requirements  of  a  new  computer  system  (GE  635),  several  large 
entities  within  the  main  program  were  separated  out  as  subroutines. 

(2)  A  minor  refinement  was  made  to  the  calorimeter  integral  calculation  to  account  for 
the  slight  variation  of  effective  view  factor  with  angle  in  a  plane  normal  to  the  arc 
centerline. 
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(b)  side  view,  section  (c)  end  view 


Figure  1.  Radiant  Interchange  Model 


(3)  A  capability  for  computing  spectral  line  emission  was  added.  , 

B.  SPECTRAL  LINES 

At  the  very  high  pressures  which  are  of  interest  to  us,  the  spectral  lines  are  very  much 
broadened  and  all  but  a  few  outstanding  ones  are  submerged  in  the  continuum.  Some  remain 
significant  however,  and  to  determine  just  how  significant,  a  few  which  can  be  clearly  distin¬ 
guished  experimentally  and  which  are  well  known  theoretically  are  chosen  for  exploration  with 
program  ARCRAD  nib.  Those  chosen  are  shown  in  Table  I;  they  include  four  hydrogen  and 
two  oxygen  lines,  the  latter  two  being  of  just  about  the  same  importance  whether  the  plasma 
is  air  or  Delrin. 

Spectral  lines  were  accounted  for  as  an  addendum  to  the  continuum  calculation  by  taking, 
for  each  line,  a  band  of  wave  numbers  within  which  radiation  from  the  line  is  significant, 
and  subdividing  that  band*,  centered  about  the  nominal  wave  number  of  the  line.  Line 
broadening  and  shift  are  both  temperature  and  density  dependent  so  they  must  be  computed 
separately  for  each  annulus.  Within  each  wave  number  subinterval,  and  for  each  annulus, 
the  absorption  coefficient  including  that  due  to  the  continuum  for  the  line  was  calculated  at  the 
center  of  each  wave  number  subinterval  in  the  following  way.  Following  Penner^),  inte¬ 
grated  absorption  is  given  by 

7/N,\  _v  _2 

kdn  =  2.375  x  10  I— )  (1  -  e  )  f^u  [cm  ]  (1) 

-OO  '  ' 

where  is  the  number  of  atoms  in  the  lower  state  of  the  transition,  Lo,  is  the  Loschmidt 
number,  f^u  is  the  oscillator  strength,  and  the  exponential  accounts  for  stimulated  emission. 
The  assumption  of  LTE  permits  calculation  of  Njt  from  the  Boltzmann  distribution,  approxi¬ 
mated  as 


N 


-v 


Z 


Z  -v. 

L  -  -  1 


i=l 


gie 


(2) 


where  g  is  statistical  weight  and  v  =  hl^/k^T.  The  natural  line  prefile  was  used  so  that  the 
normalized  probability  distribution  for  emission  in  the  wave  number  interval  dn  is  given  by 


*  The  total  number  of  subintervals  cannot  exceed  the  71  intervals  used  for  the  continuum 
calculation. 
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(3) 
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The  half-half-width  (An)jy2  *s  ^e  wave  number  interval  within  which  the  line  intensity 
decreases  from  its  maximum  value  to  one  half  its  maximum  value.  It  is  tabulated,  along  with 
the  line  shift,  n^,  by  Griem^11^  for  many  lines,  including  the  oxygen  lines  in  Table  I.  For 
hydrogen  lines,  (An)i/2  was  determined  from  published  curves  of  hydrogen  line  intensity^). 

The  absorption  coefficient  due  to  continuum  radiation  (assumed  constant  and  equal  to  its 
value  at  the  nominal  line  center  wave  number)  must  also  be  added,  so  that  we  have,  finally 


k  (n,  T) 


k  dn]  J(n-n)  +  k  ,  (T) 
o'  cont  ' 


(4) 


Having  set  up  the  absorption  coefficients  and  wave  number  intervals,  the  calculation  of 
radiation  including  a  spectral  line  then  proceeds  in  precisely  the  same  fashion  as  the  continu¬ 
um  previously  computed.  Any  combination  of  the  lines  listed  in  Table  I  can  be  included  in 
a  computer  run. 

C.  TWO  ZONE  MODEL 

The  lower  temperature,  outer  part  of  the  arc,  contributes  relatively  little  either  to 
power  dissipation  or  to  radiant  emission  and  the  accurate  determination  of  temperature 
profiles  in  the  outer  region  is  very  difficult.  For  these  reasons,  the  arc  is  divided  into  two 
zones,  Interior  (Zone  I)  and  exterior  (Zone  II).  In  Zone  I,  once  the  temperature  profile  has 
been  established,  as  discussed  below,  the  annulus  boundaries  are  so  determined  that  the 
weighted  average  temperature  of  each  annulus  is  a  multiple  of  1000°K.  This  procedure 
results  in  a  representative  but  tractable  number  of  annuli  and  eliminates  the  need  for  temper¬ 
ature  interpolation  of  absorption  coefficient  data.  In  Zone  II,  the  profile  was  assumed  linear, 
in  1000°K  increments  down  to  1000°K.  In  the  outer  layers,  radiant  emission  is  very  small 
but  absorption  is  strong  enough  to  prevent  vacuum  UV  radiation  from  reaching  the  wall. 

Program  ARCRAD  mb  was  set  up  to  use  a  temperature  profile  of  the  form 

<T  -  TA)/(TCL  -  TA>  =  1  -  (r/rA>X  <5> 

because  the  fullness  of  the  profile  is  then  governed  by  a  single  parameter  x,  and  the  symmetry 
condition  at  the  centerline  is  satisfied  for  any  x  >  1.  Best  fit  to  a  particular  temperature 
profile  (in  the  least  squares  sense)  is  obtained  if  Tql  and  T^  are  also  optimized.  The 
standard  least  squares  procedure^®)  leads  to  three  simultaneous  non-linear  equations  in  x, 
Tql  and  T^  which  are  difficult  to  solve.  However,  at  fixed  there  are  only  two  equations, 
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Table  II  Spectral  Lines  Included  in  ARCRAD  nib 


Line 

B 

n 

o 

Half-half  Width 
(Refs.  11  and  13) 

Lower  Energy  Level  • 
(Ref.  14) 

Oscillator  Strength 
(Ref.  15) 

(h 

(cm"1) 

A 

(cm  ) 

— 

10, 000°K 

20,  000°K 

40, 000°K 

L 

a 

1215 

82259 

0.  00676 

0. 00676 

0.  00676 

0 

0.4162 

1025 

97492 

0.  0121 

0.  0121 

0.  0121 

0 

0. 0791 

Hc 

6562 

15237 

0.871 

0.  871 

0.871 

82,259 

0.6407 

H0 

4861 

20570 

4.76 

4.76 

4.76 

82,259 

0.1193 

°777 

7773 

12864 

0.  0327 

0.  0443 

0.566 

73,768 

0. 922 

°845 

8445 

11839 

0.0528 

0.  0708 

0. 0895 

76,795 

0.898 

both  of  which  can  be  solved  for  dimensionless  TqL  as  a  function  of  x.  If  6.  =  (T.  -  TA)/TA, 
and  the  temperature  is  tabulated  for  m  values  of  radius,  these  equations  are 
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CL 


m 

s  9.  [l  -  (r./r/j 
i=l 

m  2 

E  [l-(r  /r  )X] 
i=l  1  A 


0 


CL 


m 

S,  ®1  (rl/rA)X  in  <rl/rA> 
_ 1=1 _ 

m 

2  (r/rA)X  ( in  r/rA)  Cl-  (r/rA)X  3 

i=l 


(6) 


(7) 


A  computer  program  providing  automatic  successive  trial  values  and  linear  interpolation  in 
the  neighborhood  of  the  solution  led  to  rapid  convergence.  Repetition  for  several  values  of 
Ta  gave  the  required  fit. 

D.  TEMPERATURE  PROFILES  FROM  NON-OPTICALLY  THIN  ARCS 

Determination  of  arc  temperature  from  measurements  of  continuum  intensity  depends  on 
an  Abel  inversion  (Appendix  I)  of  a  series  of  intensity  measurements  along  a  line  of  sight  at 
successive  increments  of  displacement,  x,  from  the  cylindrical  axis  of  symmetry.  When  the 
plasma  is  optically  thin  (negligible  self-absorption),  the  conventional  Abel  inversion  converts 
intensities,  I  (x)  to  local  emission  coefficients  i  (r),  which  can  then  be  related  to  temperatures, 
assuming  local  thermodynamic  equilibrium,  if  continuum  absorption  coefficients  are  avail¬ 
able  (see  Section  m). 


Note  that  and  c 2  are  the  Planck  radiation  constants  and,  again,  v  =  hlVkgT. 

The  method  described  herein  is  based  on  use  of  the  Abel  inversion  integral  in  an  itera¬ 
tive  way  suggested  by  Griem^11^.  Successive  approximations  to  I*(x),  the  intensity  which 
would  have  been  observed  if  no  absorption  had  taken  place,  are  computed  with  measured 
intensities  used  as  an  initial  approximation  to  I*(x).  Inversion  yields  emission  coefficients 
and  hence  a  temperature  distribution  which  can  be  used  to  determine  local  absorption.  Abel 
integration,  Appendix  I,  can  then  be  applied  to  compute  approximations  to  both  I*(x)  and  I(x). 
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If  I(x)  has  been  computed  from  the  correct  I*(x),  it  will  match  the  observed  I(x).  In  general, 
the  observed  and  computed  values  will  not  match,  so  a  new  approximation  to  I*(x)  is  obtained 
by  multiplying  each  current  local  I*(x)  by  the  ratio  of  observed  to  computed  local  I(x). 


Optical  depth  is  found  by  an  Abel  integration  along  the  line  of  sight  of  the  form  of  Equation 
6 ,  Appendix  I. 


r  (x) 


i (r) r  dr 


(9) 


The  Planck  function,  Iq<,  is  known  because  temperature  has  been  determined  from  i  (r). 
Because  i  (r)  and  I*j«  are  known  only  at  discrete  points,  the  integration  must  be  numerical,  as 
discussed  in  detail  in  Appendix  I. 

Creme rs'(^)  procedure  for  performing  the  Abel  inversion  was  followed.  The  plasma 
radius  is  divided  into  zones  of  equal  size  and  a  polynomial  is  fitted  to  the  data  in  each  zone 
by  the  method  of  least  squares^16'.  Details  of  this  process  are  also  given  in  Appendix  I  and 
the  computer  programs  developed  to  carry  it  out  are  discussed  in  Appendix  II. 
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SECTION  in 


COMPUTATION  OF/WITH  DELRIN  PROPERTIES 

To  study  electric  arc  behavior  one.  must  first  have  an  arc  that  is  well  defined  and  will 
hold  steady  long  enough  to  be  photographed;  in  addition,  the  analysis  requires  cylindrical 
symmetry.  The  practical  aspects  of  achieving  this  experimentally  are  discussed  in  Section 
V,  but,  briefly,  the  most  effective  method  we  have  found  is  stabilization  by  ablation  from  a 
constrictor  tube  made  of  an  organic  plastic.  The  obvious  drawback  is  that  we  are  then  limited 
to  a  plasma  composed  of  the  ablation  products,  but  this  is  not  serious  if  the  ablation  is  clean 
and  uniform,  and  if  the  plasma  constituents  are  known. 

An  ideal  material  for  the  purpose  is  Delrin,  a  material  which  is  pure  polyformaldehye 
(H2CO)n,  except  for  less  than  0.  5  percent  anti-oxidants  containing  C,  H,  N  and  O;  according 
to  the  Philadelphia  Product  Information  office  of  the  DuPont  Chemical  Company.  With  the 
basic  chemical  composition  known,  it  was  then  possible  to  develop  tables  of  plasma  proper¬ 
ties  for  use  in  subsequent  analyses.  Since  the  plasma  is  highly  ionized,  and  the  phenomena 
of  interest  are  closely  associated  with  free  electrons,  no  qualitative  and  little  quantitative 
difference  is  expected  from  behavior  of  arc  plasma. 

A.  THERMODYNAMIC  PROPERTIES  OF  DELRIN 

Thermodynamic  properties  were  computed  by  Browne  using  a  program 
already  developed  and  previously  applied  to  air  and  other  plasmas.  The  remainder  of  this 
subsection  quotes  from  a  report  (18)  written  about  the  air  plasma  work: 

’The  procedure  which  is  followed  ....  is  to  generate  a  set  of  ideal  gas 
thermal  functions,  that  is,  enthalpy,  free  energy  and  specific  heat,  for  each  of 
the  species  of  interest  in  the  system.  The  list  of  constituents  treated  in  these 
studies  include  the  various  atoms,  atomic  ions,  diatoms,  diatomic  ions,  poly¬ 
atomic  molecules  and  the  electron.  In  addition  to  the  thermal  functions  of  each 
specie,  the  heat  of  formation  is  required.  The  heat  of  formation  data  is  often 
inferred  from  thermochemical,  spectroscopic  or  ionization  potential  measure¬ 
ments.  The  equilibrium  composition  of  the  mixture  is  found  by  effectively 
minimizing  the  Gibbs  free  energy  given  any  two  state  variables  which  describe 
the  thermodynamic  state  of  the  system.  The  mixture  calculations  are  predicted 
on  the  absence  of  forces  between  charged  species.  Kncwling  the  equilibrium 
composition,  the  thermodynamic  properties  of  the  mixture,  namely,  enthalpy, 
entropy,  etc. ,  can  then  be  readily  determined . 

'The  calculation  procedure  ....  involved  the  evaluation  of  the  translational 
and  electron  partition  functions  and  their  contribution  to  the  thermal  functions. 

The  electronic  energy  level  summation  was  extended  over  all  energy  levels  up  to  and 
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(14) 

beyond  the  ionization  level  as  tabled  in  Moore  '  '  .  This  procedure  has  been 
demonstrated  to  yield  essentially  the  same  thermodynamic  properties  for  atoms 
and  atomic  ions  at  temperatures  below  15,  OOO^K,  as  predicted  by  using  one  of 
the  more  exact  electronic  cutoff  procedures  '  . 

"The  ideal  gas  thermal  functions  and  thermodynamic  properties  of  diatoms  and 
diatomic  ions  constitute  a  composite  of  low  and  high  temperature  calculations. 

The  low  temperature  calculations  were  based  on  the  rigid  rotator-harmonic 
oscillator  approximations  with  centrifugal  stretching  and  vibrational  anharmonicity 
corrections  included.  The  high  temperature  calculations  are  based  on  the  virial 
methodology.  Herethe  classicial  second  virial  coefficient  is  evaluated  using  the 
Morse  potential.  The  second  virial  coefficient,  in  conjunction  with  its  first  and 
second  temperature  derivative  is,  in  turn,  related  to  the  thermal  functions  of 
the  diatom  through  the  partition  function.  The  polyatomic  molecular  thermo¬ 
dynamic  properties  are  based  on  the  rigid  rotator-harmonic  oscillator  approxi¬ 
mation  . 

"The  mixture  equilibrium  composition  calculation  procedure  is  the  well- 
known  Brinkley  method  '  ’  in  which  the  mixture  is  considered  to  be  composed 
of  ideal  gases  without  charge  interaction.  The  system  of  algebraic  equations 
expressing  the  free  energy  constraints  is  solved  along  with  the  atomic  conserva¬ 
tion  statements  by  use  of  the  Newton-Rhapson  technique  .  .  .  .  " 

Results  of  the  thermodynamic  equilibrium  calculations  are  presented  as  thermodynamic 
state,  Figure  2,  and  composition  as  a  function  of  temperature  and  pressure,  Figures  3a-3d, 
Individual  species  concentrations  were  also  available  for  use  in  the  determination  of  radiation 
absorption  coefficients. 

B.  TRANSPORT  PROPERTIES 

An  examination  of  published  values  of  transport  properties  for  several  types  of  plasmas 
indicated  that  results  for  several  different  plasmas  above  15,000°K  by  one  author  tended  to 
differ  about  as  much  or  sometimes  less  than  results  for  the  same  plasma  from  different 
authors  (21-23).  The  correlations  of  Weber  (^4)  which  have  previously  been  applied  to  the 
air  plasma  were  therefore  also  used  for  Delrin  in  ARC  RAD  Illb. 

C.  RADIATION  DATA 

In  computing  absorption  coefficients  for  continuum  radiation,  radiative  processes  taken 
into  account  were:  free-bound  deionization  of  singly  ionized  H,  C,  and  O;  photodetachment  of 
electrons  from  the  negative  ions  H  ,  C  ,  and  O  ,  and  free-free  radiation  (Bremstrahlung). 
These  were  considered  to  be  the  primary  processes  since  peak  arc  temperatures  are  of  the 
order  of  20, 000°K  where  these  species  predominate.  Below  10, 000°K  an  enormously  complex 
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system  of  molecular  band  radiation  exists*  but,  as  long  as  the  plasma  layer  at  these  tempera¬ 
tures  is  not  optically  thick,  the  radiant  intensity  effect  is  negligible.  For  this  reason  only  the 
C>2  Schumann-Runge  system  was  included  as  a  sample.  As  is  the  case  for  air,  molecular 
absorption  of  vacuum  ultraviolet  radiation  is  an  important  process  at  low  temperatures  since 
this  serves  to  prevent  such  radiation  from  reaching  the  arc  column  boundary. 

Sources  and  methods  are  summarized  below. 


a.  H  Free  Bound 

For  the  hydrogen  atom  a  direct  calculation  of  the  deionization  continuum  is  possible  2(^ 
Expressed  as  an  absorption  cross-section  the  equation  is 


cth  - 


2  6  3 

32  n  e  Ry  exp  (-x, )  K 

r  '  1  mks 


3y/3  h3  c3  n3 


10  .  . 
exp  (x.) 

** *  O 

n.<n  i 
*-  l 


exp  (x  -  1) 

—  +  - - - 

2x, 


2 

(cm  ) 


(10) 


where  Ry  is  the  Rydberg  constant  and 


x.  = 


_  h  c  Ry 
kB 


(11) 


x. 

l 


.2 


(12) 


n.  = 


SSL 

2 

l 


(13) 


b.  O  and  C  free  bound 

The  quantum-defect  method  of  Burgess  and  Seaton  ^27^  was  used  by  Sherman  and  Kulander 
to  compute  cross  sections  for  radiative  recombination  (absorption  coefficient  per  particle,  in 
cm2)  of  several  species  including  O  ^28^  and  C  Typical  cross  sections  are  shown  in 

Figure  4.  Although  values  were  available  only  to  20,000°  K,  they  change  very  slowly 

(25) 

*  Main  and  Bauer  '  list  44  different  systems  involving  C,  H,  and  O. 
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with  temperature,  so  a  second-order  LaGrange  interpolation  was  used  to  extrapolate  to 
25, 000°K  based  on  values  at  16,000°,  18,000°,  and  20,000°K. 

c.  H- 


Radiation  associated  with  the  negative  hydrogen  ion  has  been  the  subject  of  considerable 
interest  to  astrophysicists  for  many  years.  Chandresekhar  (30>  31 )  gives  a  curve  (Figure  5a) 
for  the  cross  section  for  deionization  and  a  tabulation  of  both  free-free  and  free-bound  radia¬ 
tion  as  a  function  of  temperature.  His  table  of  the  absorption  coefficient  per  neutral  hydrogen 
atom  and  per  unit  electron  pressure  for  free-free  radiation  associated  with  the  negative 
hydrogen  ion  was  found  to  be  well  represented  by  a  double  logarithmic  interpolation  of  selected 
values  (Table  II). 


Table  II.  Selected  Values  of  k 


Wave  Number 
(cm-1) 

0  =  5040/T  (°K) 

2 

0.5 

2000 

8.4  x  10-25 

3.  6  x  10-25 

8000 

4. 12  x  10-26 

8.  9  x  10-27 

20000 

8.  9  x  10~27 

1.27  x  10-27 

The  quantity  k  is  defined  as  absorption  coefficient  per  neutral  hydrogen  atom  and  per 
unit  electron  pressure  (atm),  (1500  cm-1  <  n  <  25, 000  cm-1).  Linear  absorption  coefficient 
is  then  given  by 


k<V 


ill)  (\)  (to) 


x  1.  013  x  10^  (cm  ^) 


(14) 


d.  C 

Cross  sections  have  been  calculated  both  by  Breene  and  by  Kulander  and  Sherman  (28\ 
The  latter’s  values  were  used.  They  are  based  on  the  quantum  defect  method  and  are  more 
than  an  order  of  magnitude  larger  (Figure  5b). 


e.  O 


Cross  sections  have  been  measured  by  Branscome  et  al  (33)  up  to  about  3  ev  (24, 000 
cm-1).  Above  that  wave  number  the  procedure  of  Meyerott  as  discussed  in  Penner  (12'  was 
followed,  with  the  constants  adjusted  to  fit  Branscome’s  more  recent  measurements  (Figure  5a). 
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f.  Bremsstrahlung 


The  Kramers-  Unsold  equation  was  used,  with  Gaunt  factor  and  ion  charge  both  at 
unity(34). 


kff 


3,708  x  10 
3  3 


8 


n 3  yfF 


(1  -  e~V)  N  ( E  N  )  i  =  C+,  H+,  0+ 

l 


(15) 


g.  02  Schumann-Runge 

q  o 

Data  were  taken  from  curves  for  spectral  radiance  J(watt/cm  -str-A)  published  by  Breen 
and  Nardone(^)  and  converted  to  absorption  coefficients. 

h.  Vacuum  UV  Molecular  Absorption 

The  predominant  molecular  species  in  Delrin  plasma  are  HgO,  CO,  CO„,  and  CH^. 
Absorption  by  these  species  was  accounted  for  by  averaging  the  experimental  data  collected 
by  Schultze,  et  al(36). 


Line  radiation  absorption  coefficients  are  discussed  in  Section  II  part  B,  with  key  data 
shown  in  Table  I. 


Figures  6a-6e  show  continuum  spectral  radiance  of  Delrin  at  atmospheric  density  and  at 
5000  K  temperature  increments.  The  graphs  are  computer  plotted  for  a  1  millimeter  thick 
slab  and  show  the  individual  species,  total  spectral  radiance  and  the  black  body  limiting  curve. 
As  expected,  the  Bremsstrahlung  and  deionization  continua  are  dominant  and  the  vacuum 
ultraviolet  region  is  nearly  black  body  even  at  low  temperature. 

D.  COMPUTATION  SUMMARY 

The  block  diagram,  Figure  7,  summarizes  the  computations  performed  on  this  project. 
The  absorption  coefficients  and  thermodynamic  data  already  discussed  were  inputs  to  program 
ARCRAD  mb.  They  were  also  the  source  for  calculations  of  emission  coefficient  as  a  function 
of  temperature  at  a  specified  wave  number,  to  be  used  in  the  temperature  profile  measurements 
discussed  in  Section  HD  and  part  C  of  Appendix  II. 

Figure  8  shows  typical  emission  coefficient  behavior.  At  the  infrared  wavelength  chosen, 
Bremsstrahlung  and  radiation  are  dominant,  as  can  be  seen  from  the  dashed  curve  on  the 
figure  which  was  calculated  directly  from  Equations  10  and  15. 
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Figure  3b.  Delrin  Composition,  P  =  30  atm 
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Figure  3c.  Delrin  Composition,  P  =  100  atm 
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Figure  3d.  Delrin  Composition,  P  =  200  atm 
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ABSORPTION  CROSS  SECTION  (CM*) 


2 


ABSORPTION  CROSS  SECTION  (CM  ) 


Figure  4b.  Carbon  Deionization  Absorption  Cross  Section 
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Figure  5a.  Photoinization  Cross-Section,  H  and  O 
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Figure  5b.  Photoinization  Cross  Section,  C  ,  Two  Sources 
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Figure  6a.  Delrin  Spectral  Radiance,  1  mm  Slab,  T  =  5000°K, 
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Figure  6b.  Delrin  Continuum  Spectral  Radiance,  1  mm  Slab,  T  =  10,000°D,  p/pQ  =  1 
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Figure  6c.  Delrin  Continuum  Spectral  Radiance,  1  mm  Slab,  T  =  15,000°K,  p/p  =  l 
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Figure  6d.  Delrin  Continuum  Spectral  Radiance,  1  mm  Slab,  T  =  20,000  K  p/p 
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Figure  6e.  Delrin  Continuum  Spectral  Radiance,  1  mm  Slab,  T  =  25, 000°K  p/p 


Figure  7.  Block  Diagram  of  Computations 


EMISSION  COEFFICIENT  (lO9  ERG/SEC“CM3-&-STR) 


Figure  8.  Emission  Coefficient  for  use  in  Temperature  Profile  Determination 
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SECTION  IV 


APPARATUS 

The  conceptual  design  of  the  experimental  apparatus  remained  largely  unchanged  from 
last  year.  A  brief  description  will  follow,  but  for  further  details  see  Reference  1.  Figure  9 
shows  an  overall  view  of  the  experimental  area  which  is  enclosed  in  armor  plate  for  safety 
protection.  The  ensemble  has  been  rearranged  since  last  year. 

Arcs  were  maintained  for  only  a  few  milliseconds  because  the  large  heat  flux,  of  the 
order  of  100  kw  /Cm2,  inevitably  produces  wall  ablation  after  a  millisecond.  Even  with  an 

ablation  type  constrictor,  it  is  undesirable  to  run  for  times  longer  than  about  10  milliseconds 
because  constrictor  diameter  increases  rapidly  and  this  constitutes  a  time  varying  boundary 
condition  with  strong  effect  on  arc  parameters.  Power  is  supplied  by  shaping  the  discharge 
of  a  500  /if,  18,  000  volt  capacitor  bank  to  square  pulse.  Capability  is  250  amperes  for 
10  milliseconds  with  correspondingly  higher  currents  for  shorter  times.  Switching  is  by 
triggered  ball  gaps. 

3 

The  aluminum  test  chamber  is  rated  at  6,  000  psi  and  has  a  volume  of  220  cm  .  Optical 
observation  is  made  through  two  symmetrically  placed  Plexiglas  windows,  1.4  x  4  cm. 

Synchronization  of  the  arc  and  spectrographic  exposure  or  other  measurements  is 
achieved  with  shuttered  light  sensors  positioned  around  a  rotating  mirror.  The  same  mirror 
serves  as  the  central  time  resolving  element  of  the  spectrograph  optical  system.  An  exposure 
of  100  microseconds  was  a  successful  compromise  of  the  conflicting  requirements  of  plate 
exposure  and  stopping  the  remaining  motion  of  the  arc. 

Spectra  were  recorded  on  Kodak  IN  and  IF  plates  with  a  Hilger  medium  glass  spectro¬ 
graph.  Photographic  calibration  was  made  with  the  anode  of  a  carbon  arc  as  intensity 
standard  (37,  38).  Identical  optical  system  and  exposure  time  were  used  for  both  arc  and 
intensity  standard. 

Electrical  measurements  were  made  with  CRT  oscilloscopes.  A  Tektronix  502A  was 
used  for  recording  the  weak  signal  of  the  calorimeter  thermocouple  and  Tektronix  551  and  545 
units  for  all  other  recordings. 

Measurements  were  taken  at  100,  150,  and  200  atmospheres  initial  pressure.  The 
pressure  rise  in  the  test  chamber  during  arcing  was  measured  with  a  Kistler  #602A  quartz 
pressure  transducer  so  that  the  pressure  at  the  time  of  the  spectral  exposure  was  known 
exactly. 
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Figure  9a.  Overall  view  of  experimental  ensemble.  The  arc  is  housed 
in  aluminum  pressure  chamber  at  center.  Two  windows  and  a  beam  splitter 
allow  simultaneous  observation  of  the  arc  by  spectograph  and  framing 
camera  at  right  and  streak  camera  at  left. 


Figure  9b.  Close-up  view  of  spectrograph  optics.  A  magnified  image  of  the 
arc  is  focused  on  rotating  mirror  at  left;  a  second  image  is  formed  on  the 
spectrograph  slit  behind  the  shutter.  The  boxes  grouped  around  the  mirror 
contain  a  light  source,  light  sensors,  and  electronics  to  synchronize  the 
operation  of  the  rotating  mirror,  the  spectrograph  shutter,  and  the  arc.  The 
flip-flop  mirror  (partially  hidden  by  arc  chamber)  is  positioned  midway 
between  the  arc  and  the  carbon  arc  intensity  standard  so  that  both  may  be 
photographed  by  the  same  system. 


SECTION  V 


ARC  CONSTRICTORS  AND  ARC  STABILITY 

Pertiaps  the  most  important  and  critical  part  of  the  experimental  apparatus  is  the  arc 
constrictor.  Our  problems  with  constrictor  design  and  some  important  observations  shall 
therefore  be  described  in  the  following  separate  chapter. 

A.  CONSTRICTOR  STABILIZATION 

Previously,  two  types  of  constrictors  were  used.  ^  One  was  the  Maecker  type  constrictor, 
which  consists  of  a  stack  of  copper  discs,  forming  a  cylindrical  tube  for  wall  stabilization. 

The  other  was  an  ablation  type  constrictor  with  a  cylindrical  wall  made  of  Plexiglas  which 
vaporized  under  action  of  the  arc,  generating  a  cylindrically  symmetric  flow  field  of  wall 
material. 

High  speed  photographic  work  with  Dynafax  framing  and  streak  cameras  had  shown  that 
the  Maecker  constrictor  could  not  generate  stable  arcs  at  high  specific  heat  fluxes  but 
exhibited  strong  kink  instabilities,  with  the  instability  increasing  with  increasing  pressure 
and/or  arc  energy.  Whereas  it  was  still  possible  to  measure  the  radial  temperature  distri¬ 
bution  of  a  100  atmospheres,  100  amperes  are  with  an  arc  energy  of  20  kw/cm,  an  arc  at 
100  atmospheres,  with  400  amperes  and  80  kw/cm,  was  too  unstable  for  a  measurement 
of  the  radial  temperature  distribution. 

In  contrast  to  this,  the  ablation  type  arc  had,  in  some  preliminary  experiments,  shown 
excellent  stability.  There  is,  of  course,  the  drawback  that  the  ablation  type  arc  cannot  be 
burned  in  an  atmosphere  of  the  pressurizing  gas,  but  only  in  the  evaporation  products  of  the 
tube  wall.  To  get  N2,  or  C02,  plasmas,  etc. ,  with  an  ablation  type  constrictor,  one  must 
therefore  use  these  gases  in  solidified  form.  This  has  been  tried  in  the  case  of  C02(6);  but 
the  required  technique  is  not  simple. 

In  view  of  these  difficulties,  an  additional  effort  was  made  to  obtain  stable  air  arcs  in 
the  Maecker  type  constrictor  100  atmospheres  and  200  amperes.  This  attempt  was  not 
successful.  The  use  axial  magnetic  fields  up  to  2,000  gauss  to  enforce  stability  was  also 
to  no  avail. 

B.  RADIATION  STABILIZATION 

At  the  suggestion  of  Dr.  H.  Hurwitz,  *  the  arc  radiation  rather  than  tube  walls  was 
examined  for  defining  a  stable  channel.  A  cylindrical  mirror  around  the  arc  axis  would 
redirect  arc  radiation  into  the  arc  channel  proper  and  be  partially  reabsorbed  there. 


*  Dr.  Hurwitz  is  with  the  General  Electric  Company  Research  and  Development  Center, 
Schenectady,  N.  Y. 
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This  process  would  constitute  a  feedback  between  the  arc  and  its  radiation  and  would  help  to 
burn  the  arc  on  the  axis  of  the  cylindrical  mirror.  The  idea  was  tried  and  generated  a  plasma 
qualitatively  much  different  than  obtained  in  copper  constrictors.  The  plasma  did  not  constrict 
in  a  narrow  channel,  but  instead  continuously  expanded  during  the  1  millisecond  duration  of 
the  experiment,  Figure  10.  While  the  idea  of  a  radiation  stabilized  arc  shows  promise  and 
may  well  have  important  applications,  the  qualitatively  very  different  character  of  the  dis¬ 
charge  would  require  a  major  amount  of  work  before  it  could  be  used  for  measurements  in 
the  present  study.  Therefore,  we  did  not  follow  this  interesting  idea  in  more  detail,  but 
turned  with  a  major  effort  to  the  ablation  type  constrictor  which  had  previously  shown  good 
stability. 

C.  DEVELOPMENT  OF  THE  ABLATION  TYPE  CONSTRICTOR  FOR  OPTICAL  MEASUREMENTS 

The  stability  of  the  ablation  type  constrictor  turned  out  be  be  excellent — as  long  as  there 
were  no  openings  in  the  walls  of  the  cylindrical  tube.  Figure  11  shows  some  high  speed  camera 
frames  and  the  simultaneous  streak  picture  of  an  arc  in  a  Plexiglas  tube.  After  a  very  short 
unstable  period,  the  arc  steadies  and  remains  steady  during  the  complete  duration  of  the 
experiment.  However,  the  simple  cylindrical  tube  of  Plexiglas  has  some  important  draw¬ 
backs  if  side-on  temperature  measurements  are  to  be  made:  (1)  the  tube  constitutes  a  bad 
cylindrical  lens  which  distorts  the  arc  image;  (2)  Plexiglas  evaporates  only  partially  under 
action  of  the  arc  and  shows  some  surface  carbonization  (The  carbon  particles  often  are 
focussed  on  the  spectrograph  slit  and  cause  a  streaky  spectrogram);  (3)  The  Plexiglas  wall 
does  not  evaporate  completely  homogeneously  but  attains,  after  a  while,  a  frosty  glass 
appearance.  None  of  these  effects  are  acceptable  and,  therefore,  a  good  window  in  the 
Plexiglas  wall  was  needed  for  better  quality  optical  observation. 

The  simplest  window  would  be  a  slot  between  two  short  Plexiglas  tubes.  Figure  12. 

With  this  design,  part  of  the  plasma  streams  through  the  slot  and  the  luminous  volume  is 
larger  than  the  arc  channel  proper.  Moreover  we  found  arcs  burning  violently  unstable  in 
this  design,  apparently  because  small  maladjustments  between  the  two  Plexiglas  tubes  caused 
an  unsymmetric  and  fluctuating  flow  field. 

After  several  intermediate  steps,  we  came  finally  to  the  design  shown  in  Figure  13. 

This  design  has  a  flat  quartz  window,  relatively  far  removed  from  the  arc  channel.  The 
window  is  mounted  tightly  to  the  body  of  the  plastic  "tube"  to  prevent  plasma  from  streaming 
out  into  the  window  area.  The  design  is  symmetric  with  two  windows  at  180°  for  simultaneous 
optical  observation  with  our  time  resolved  spectrograph  and  a  streak  camera.  As  "tube" 
material,  Delrin  was  chosen,  both  for  its  simple  composition  and  for  its  "non-carbonization" 
characteristic.  The  flat  window  was  Suprasil  quartz,  chosen  for  minimum  energy  absorbtion 
of  the  extremely  strong  far  ultraviolet  radiation  of  our  arc. 

This  constrictor,  which  will  be  referred  to  as  our"optical"  constrictor,  maintains  a 
relatively  stable  arc,  as  shown  by  Figure  14.  The  figure  reveals,  however,  that  the  arc 
does  not  burn  strictly  on  the  axis  of  the  tube  but  slightly  to  one  side,  and  some  very  slow,  low 
amplitude  motions  persist.  This  degree  of  stability,  however,  is  sufficient  for  our  very 
short  time  spectrographic  exposures. 
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Besides  the  "optical"  ablation  type  constrictor,  several  others  were  built  for  measure¬ 
ments  of  electrical  gradients,  plasma  pressure,  and  radial  heat  flux.  For  these  measure¬ 
ments,  the  rather  complicated  optical  window  was  not  necessary  and  therefore  omitted. 

All  constrictors,  however,  were  basically  Delrin  tubes  with  minor  design  modifications  as 
demanded  by  the  specific  measurement.  They  are  described  in  the  sections  where  we 
report  the  respective  measurements. 
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Figure  10  b.  Magnified  views  of  the  "radiation  stabilized"  arc  suggested  by  Hurwitz.  The 
port  is  10mm  x  5mm.  The  300  ampere  arc  was  initiated  by  exploding  a  0.  00025  inch  diameter 
tungsten  fuse  wire  which  was  accurately  centered  at  the  focus  of  the  polished  stainless  steel 
tube.  Compare  this  diffuse  arc  with  the  constricted  ones  elsewhere  in  this  report  (Figure  11) 
and  last  year  (Ref.  1,  Figure  36). 
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Figure  11.  Stability  of  a  Plexiglas  ablation  arc  as  seen  by  a  streak  camera  and  a 
framing  camera.  Note  that  the  arc  is  initially  unstable  but  becomes  stable  in  a 
few  hundred  microseconds.  The  arc  in  the  streak  camera  film  is  magnified  5.5  x; 
in  the  frame  pictures  magnification  is  2x.  The  omitted  frames  are  completely 
stable  and  identical  with  frames  #3,  10  and  20. 


Figure  12.  Schematic  design  of  a  simple,  two-piece  ablation  constrictor  with  provision 
for  spectroscopic  observation  of  the  arc.  Such  a  device  was  found  to  produce  an  unstable 
arc  because  of  slight  misalignment  of  the  two  pieces  and  enlargement  of  the  arc  in  the 
free  region.  For  final  design  see  Figure  13. 
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Figure  13  a.  Schematic  design  of  optical  constrictor.  The  arc  tube  is  made  from  one  piece 
with  minimal  openings  for  observation.  The  Suprasil  quartz  windows  provide  a  pressure 
seal  yet  pass  almost  all  of  the  incident  radiation. 


Figure  13  b.  Photograph  of  an  assembled  optical  constrictor  and  parts.  Assembly  was 
by  glue  placed  away  from  the  arc  to  avoid  ablation. 
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Figure  14.  Streak  camera  traces  of  a  stable  and  an  unstable  test  in  Delrin  constrictors.  Minor 
design  changes  can  destabilize  the  arc  and  make  measurements  impossible.  The  view  slit  width  must 
be  the  arc  tube  radius;  the  view  slit  must  be  sealed;  and  the  "dead  volume"  between  the  window  seal 
and  the  arc  tube  must  be  less  them  or  equal  to  minimal. 
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MEASUREMENT  OF  ARC  PARAMETERS 

It  is  of  considerable  importance  to  check  whether  the  arc  is  uniform  along  its  length.  If 
the  arc  is  axially  nonuniform,  observations  would  be  required  at  a  number  of  points  along 
the  arc  axis,  which  is  time  consuming.  Plasma  flow  from  the  constrictor  ends  was  in  the 
range  of  1  to  2  grams /second,  which  might  be  expected  to  produce  a  large  pressure  rise  in 
the  constrictor  and  thus  an  undesirable  axial  variation  of  other  arc  properties.  Meaningful 
calculation  was  difficult  because  most  of  the  axial  flow  occurs  near  the  walls  where  little  is 
known  of  the  temperature  distribution,  and  there  is  also  the  possibility  that  a  thin  layer  of 
Delrin  liquid  flows  along  the  walls.  Consequently,  direct  measurements  of  the  pressure  rise 
at  the  center  of  the  constrictor,  where  spectral  observations  are  made,  was  deemed  advis¬ 
able,  along  with  voltage  gradient  measurements  for  several  constrictor  lengths. 

A.  PRESSURE  RISE  IN  THE  CONSTRICTOR 

For  the  pressure  measurement  a  special  arc  constrictor  was  built  as  shown  in  Figure  15. 
A  Delrin  port  connects  the  arc  duct  to  a  steel  tube  of  0.  009  inch  inside  diameter,  which  leads 
to  a  quartz  pressure  transducer  (Kistler  type  601  A).  This  arrangement  enables  us  to  read 
the  pressure  increase  in  the  center  of  a  constrictor  during  arc  operation.  The  system  was 
tested  for  rise  time  by  recording  the  transducer  output  while  an  over-pressure  protection 
diaphragm  on  the  test  chamber  was  intentionally  burst.  A  simultaneous  record  of  pressure 
drop  as  taken  from  a  transducer  directly  connected  to  the  pressure  chamber  was  used  as  a 
comparison  standard.  From  this  measurement  overall  system  response  time  was  found  to 
be  of  the  order  of  1  millisecond,  adequate  for  our  purposes. 

The  pressure  transducer  was  then  used  to  measure  rise  of  the  plasma  pressure  in  the 
constrictor  during  arcing.  It  was  found  that  the  rise  was  very  nearly  the  same  as  that  for 
the  entire  high  pressure  container:  approximately  three  atmospheres  during  a  5  millisecond 
run  at  100  atmospheres.  With  such  a  small  pressure  drop  in  the  tube,  there  should  be  no 
major  effect  from  the  flow  field  on  the  structure  and  energy  balance  of  the  arc. 

B.  VOLTAGE  GRADIENT  AND  AXIAL  INTENSITY  VARIATION 

Figure  16  shows  the  experimental  arrangement  used  to  determine  voltage  gradients.  The 
Delrin  constrictor  is  a  simple  cylinder  placed  between  carbon  electrodes  accurately  spaced 
1.  5  mm  from  each  end.  The  length  of  the  constrictor  was  varied  in  successive  shots  and  the 
arc  voltage,  read  from  differential  probes,  was  plotted  as  in  Figure  17  to  give  the  voltage 
gradient.  Despite  the  slight  scatter  in  the  data  points  it  is  clear  that  there  is  no  systematic 
deviation  from  linearity  in  the  range  of  lengths  tested.  The  non-zero  voltage  at  "zero" 
length  corresponds  to  the  potential  drops  in  the  spaces  between  constrictor  and  electrodes 
plus  the  anode  and  cathode  drops.  The  constrictors  used  for  spectral  investigation  were  16 
mm  long,  as  were  the  constrictors  tested  for  pressure  rise. 
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The  constancy  of  the  voltage  gradient  along  the  tube  agrees  well  with  our  earlier  obser¬ 
vation  of  a  negligible  pressure  gradient  along  the  arc  and  also  with  the  observation  that 
Plexiglas  arcs  show  no  intensity  variation  along  the  tube  length.  It  appears,  therefore,  that 
at  the  conditions  of  our  experiments,  no  significant  axial  inhomogeneity  exists  in  the  ablation 
type  arcs. 

Figures  18  and  19  show  the  results  of  gradient  measurements  at  higher  pressures  with  a 
similar  result.  Clearly  the  voltage  grows  with  pressure,  as  shown  in  Figure  20-  Figure  21 
shows  the  corresponding  curve  made  last  year  in  air  but  unreported  at  that  time.  Notice 
the  somewhat  weaker  dependence  on  air  pressure  as  compared  with  the  Delrin  results. 
Plexiglas  constrictors  gave  the  same  general  result  but  with  lower  voltage  gradients  (Figures 
22  and  23).  All  of  the  tests  discussed  above  were  made  at  250  amperes,  but  Figure  24  shows 
that  the  voltage  gradient  is  only  slightly  affected  by  a  change  in  current  near  this  value. 

C.  CONSTRICTOR  ABLATION  RATE 

The  ablation  rate  was  determined  by  weighing  the  constrictors  before  and  after  testing 
on  an  analytical  balance.  Losses  were  generally  several  milligrams  and  could  therefore  be 
determined  with  fair  accuracy.  Arc  energy  was  simply  computed  as  the  time-integrated 
product  of  arc  current  and  voltage;  anode  and  cathode  voltage  drops  were  neglected. 

It  was  found  that  the  ablation  rate  increased  when  the  arc  was  burned  for  longer  times 
(Figure  25).  A  rising  curve  for  short  times  could  be  expected  from  the  initial  period  required 
to  initiate  ablation.  Evidence  for  the  existence  of  this  "initiation  period"  comes  from  high 
speed  photographs,  (Figure  14)  which  show  the  arc  initially  unstable,  indicating  that  ablation 
stabilization  has  not  begun,  and  from  the  arc  voltage  records  which  show  an  initial  period  of 
high  voltage  that  drops  when  the  photographs  indicate  stability. 

Figure  25  also  shows  that  the  ablation  rate  increases  at  greater  current.  Pressure  is 
also  a  strong  factor  in  determining  the  ablation  rate.  Figure  26  shows  a  substantial  increase 
in  ablation  of  Delrin  as  pressure  rises  from  150  to  2900  psig. 
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Figure  15a.  Schematic  diagram  of  setup  for  measurement  of  pressure  rise  at  center  of  Delrin  ablating  constrictor. 
The  steel  tube  is  set  back  sufficiently  far  from  the  arc  for  "arc-over”  protection.  Careful  tests  were  made  to 
insure  that  the  response  time  of  the  system  was  adequate  for  accurately  measuring  the  pressure  rise  producted 
by  a  10  msec.  arc. 


Figure  15b.  Photograph  of  pressure  measuring  constrictor.  The  constrictor  is  shown  mounted 
in  a  Plexiglas  case  with  carbon  electrodes  positioned  for  firing.  The  copper  mesh  "coolers",  Fig.  16, 
have  been  removed  to  show  the  electrodes.  In  the  foreground  is  the  flexible  joint  of  the  steel 
tube  made  with  a  jacket  of  polyethylene  tubing. 
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Figure  16.  Delrin  constrictor  ready  for  firing.  The  plasma  exhaust  from  the  ends  o  the 
constrictor  passes  through  layers  of  a  copper  mesh  which  cool  the  plasma  and  reduce  the 
pressure  rise  in  the  arc  chamber.  The  spacing  between  the  ends  oi  the  constrictor  a  id  the 

electrodes  is  fixed  at  1.5  mm. 
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Figure  17.  Arc  voltage  for  Delrin  constrictors  of  different  lengths.  The  data  points  show  no 
systematic  deviations  from  linearity  in  the  range  of  lengths  tested.  The  slope  gives  a  value 
of  395  volts /cm  for  the  arc. 
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Figure  18.  Arc  voltages  for  Delrin  constrictors  of  different  lengths. 


CONSTRICTOR  LENGTH  (MM) 

Figure  19.  Arc  voltages  for  Delrin  constrictors  of  different  lengths. 
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GRADIENT  (V  /CM) 


Figure  20.  Voltage  gradient  of  Delrin  ablation  arc  vs  pressure 


PRESSURE  (PSIA) 


Figure  21.  Voltage  gradient  of  air  arc  vs  pressure 
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CONSTRICTOR  LENGTH  (MM) 


Figure  22.  Arc  voltage  for  different  lengths  of  Plexiglas  constrictors.  Compare  with 
Figures  17,  18  and  19. 


PRESSURE  (PSIA) 


Figure  23.  Voltage  gradient  of  a  Plexiglas  arc  as  a  function  of  pressure.  The  Plexiglas 
arc  has  a  potential  gradient  between  those  of  air  and  Delrin. 


50 


CURRENT  (AMPERES) 


Figure  24.  Voltage  gradient  of  a  Delrin  ablation  arc  vs  current.  Despite  the  scatter  of  data 
points,  it  can  be  clearly  seen  that  arc  current  has  small  effect  on  the  voltage  gradient  above 
100  amperes.  These  experiments  were  made  in  a  1. 5  mm  ID  tube  at  a  pressure  of  150 
atmospheres.  Earlier  experiments  with  2.  5  mm  ID  tubes  showed  voltage  gradients  25 
percent  lower. 
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ABLATION  RATE  (X  io"5  GRAMS/JOULE) 


ARCING  TIME  (MILLISECONDS) 


Figure  25.  Ablation  of  Delrin  at  100  atmospheres.  The  ablation  rate  increases  slowly  with 
arc  current  and  with  longer  burning  arcs.  These  results  were  taken  with  2. 5  mm  ID  con¬ 
strictors  which  were  weighed  before  and  after  the  test  on  an  analytical  balance. 
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PRESSURE  (PS1A) 


Figure  26.  Variation  of  ablation  rates  with  pressure.  The  Delrin  ablation  rate  is 
strongly  influenced  by  pressure  changes. 


SECTION  vn 


MEASUREMENT  OF  RADIANT  HEAT  FLUX 

The  radiant  heat  flux  calorimeter  consists  of  a  small  copper  ring  of  known  mass  con¬ 
centric  with  the  arc.  The  temperature  rise  in  the  ring  is  measured  with  a  copper-constrictor 
thermocouple.  The  ring  is  recessed  and  protected  from  conducted  heat  of  the  arc.  A  previous 
report  ^  gives  a  full  description  of  the  calorimeter  used  for  cascade  air  arcs,  but  several 
uesign  modifications  were  required  to  adapt  the  device  to  measurement  of  an  ablation  arc. 

Figure  27  shows  the  final  design  of  our  ablation  type  arc  calorimeter.  A  light  baffle 
was  added  to  limit  stray  reflections  and  more  closely  approximate  the  calculated  view  factor 
and  the  calorimeter  ring  is  protected  from  the  arc  plasma  by  a  Suprasil  quartz  tube.  The  tube 
eliminates  heating  by  conduction  and  convection  yet  nasses  almost  90  percent  of  the  incident 
arc  radiation.  Operation  of  the  arc  for  5  milliseconds  did  not  significantly  damage  the  quartz. 

A  special  black  Delrin  (doped  with  0. 3  percent  Carbon)  was  used  to  prevent  energy  trans¬ 
mission  through  the  normally  translucent  (40  percent  of  visible  light  through  a  1  mm  layer) 
tubing. 

In  spite  of  these  precautions,  no  useable  results  were  obtained  from  calorimeter  measure¬ 
ments.  Temperature  rise  was  generally  much  too  high  to  be  solely  due  to  arc  radiation  and 
changing  too  rapidly  to  be  reasonably  extrapolated  back  to  time  of  arc  cutoff.  *  Since  the 
expected  rise  was  about  10°  C  (in  3  milliseconds),  small  effects  could  distort  the  measurement. 


*  Vacuum  switches  turn  on  the  thermocouple  circuit  approximately  12  milliseconds  after 
arc  extinction  to  protect  it  from  the  very  high  voltages  during  arc  operation. 
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Figure  27b.  Photograph  of  Calorimeter  Constrictor  and  Parts 


SECTION  vni 


TEMPERATURE  MEASUREMENT.^ 

Arc  temperature  was  measured  spectroscopically  with  a  time  resolution  of  about 
100  microseconds.  The  time  resolving  element  was  a  rotating  mirror;  the  basic  spectro¬ 
graph  was  a  Hilger  medium  glass  instrument.  Details  of  this  apparatus  and  the  general 
measuring  techniques  used  are  discussed  in  the  previous  report  C1). 

In  view  of  the  fact  that  the  high  pressure  spectra  displayed  a  largely  continuous 
character  with  the  lines  strongly  broadened  and  even  overlapping,  no  great  variety  of 
temperature  measuring  methods  was  available.  The  plasma  was  strongly  reabsorbing 
even  in  the  continuum,  so  optically  thin  methods  could  not  be  used.  Line  radiation  emerging 
from  the  discharge  showed  little  contrast  to  the  strong  continuous  radiation.  Therefore, 
the  absolute  value  of  the  continuous  radiation  was  considered  to  be  the  most  likely  candidate 
for  a  temperature  measurement.  The  data  required  was  a  sequence  of  absolute  measure- 
ments-of  the  "side  on"  spectral  intensity  of  the  arc. 

o 

Preliminary  measurements  were  made  at  8330A.  However,  when  it  became  apparent 
that  line  wings  might  contribute  significantly  to  the  intensity,  all  the  plates,  including  the 
one  for  air,  were  reevaluated  at  7050$..  At  the  latter  wavelength,  the  continuum  dom¬ 
inates  in  air  and  the  only  line  of  importance  in  Delrin  is  the  Ha  line.  This  line  profile 
is  well  known,  and  the  curves  and  asymptotic  formula  of  Griem,  Kolb  and  Shen^)  were 
used  for  a  detailed  evaluation  of  its  contribution  to  total  absorption  coefficient.  As  tem¬ 
perature  increases,  the  line  tends  to  increase  both  in  intensity  and  breadth,  so  that  at 
7050$  the  line  contribution  is  negligible  at  10, 000°K,  equal  to  the  continuum  at  16,  000°K 
and  about  twice  the  continuum  at  25,000°K. 

The  print  of  a  typical  spectroscopic  plate  for  a  temperature  measurement  of  a  Delrin 
arc  of  100  atmospheres  and  250  amperes  is  shown  in  Figure  28.  To  obtain  sensitivity 
at  the  7050$  wavelength  and  also  in  the  area  of  the  strong  infrared  oxygen  lines  up  to 
9000$,  Kodak  I- N  plates  were  used.  These  plates,  however,  are  relatively  insensitive 
in  the  green,  as  seen  from  Figure  28.  Alongside  the  Delrin  arc  with  its  strong  continuum 
and  strongly  broadened  oxygen  (and  hydrogen)  lines,  the  plate  shows  the  intensity  markers 
of  the  carbon  standard^*  and  wavelength  markers  from  mercury- cadmium  and  rubidium 
lamps.  Note  that  the  well-known  hydrogen  line  is  so  strongly  broadened  that  it  can 
no  longer  be  distinguished  from  the  continuum. 

The  result  of  photometer  work  on  the  spectrum  of  Figure  28  is  shown  in  Figure  29, 
where  absolute  values  of  the  integral  spectral  intensity  I(x)  are  plotted  against  the  arc 
radius  in  millimeters.  This  curve  is  then  transformed  by  the  self-absorption  compensated 
Abelian  transformation  into  the  curve  of  the  local  spectral  intensity  (emission  coefficient)  i(r), 
Figure  30.  From  this,  together  with  the  theoretical  results  of  Section  n,  Part  C,  the 
radial  temperature  distribution  of  the  arc  is  computed  (Figure  31). 
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Similar  measurements  were  also  made  for  Delrin  arcs  of  250  amperes  at  150  and 
200  atmospheres.  The  results  of  the  temperature  measurement  at  150  atmospheres, 
where  the  self-absorption  correction  becomes  quite  large,  are  also  shown  in  Figure  31, 
but  the  200  atmosphere  arc  exhibited  too  much  stray  light  for  a  quantitative  temperature 
evaluation.  Details  of  the  numerical  evaluation,  including  the  previous  measurement  in 
air ,  are  included  as  Table  HI.  Figure  32  is  a  summary  plot  of  all  successful  profile 
measurements. 


57 


/ 


4046 


4358 


4678 

4799 

5085 

5460 

5790 

6438 


7800 

8521 


4341  H 


NOT 

y  VISIBLE 


4861  H  p 


NOT 

VISIBLE 


6562  H  a  VISIBLE 


7480 

7773 

8230 

8820 


0(  LINES 


HG-CD  CARBON  DELR1N 
AND  RB  STANDARD  ARC 


ARC  DATA 

PRESSURE 

CURRENT 

CONSTRICTOR 

DIAMETER 

TEST  NO. 


110  ATM 
250  AMPS 

1.5  MM 
DELRIN  79 


Figure  28.  Spectral  plate  for  the  measurement  of  the  temperature  of  a  high  pressure  Delrin 
arc.  Along  with  the  spectrum  of  the  Delrin  arc,  spectra  of  the  carbon  crater  (intensity 
standard)  and  of  line  radiators  (Hg-Cd  and  Rb)  for  determination  of  wavelength  are  shown. 
Spectrum  of  Delrin  arc  is  largely  continuous  with  a  few  strongly  broadened  lines  of  Oj  and  Hj 
still  marginally  recognizable. 
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EMISSION  COEFFICIENT,  t  (r)  (10B  ERG/CM3-ft-STR-SEC) 


Figure  30.  Local  emission  i(r).  Delrin  arc  shows  a  (probably  not  correct)  maximum  off 
axis.  This  maximum  is  reflected  also  in  the  radial  temperature  distribution  (Fig.  31)  but 
only  to  a  minor  degree. 
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Figure  31.  Radial  Temperature  distributions  in  High  Pressure  Delrin  Arcs  based  on  intensity  measurements  at 
7050.8. 


Table  m.  Results  of  Temperature  Profile  Measurements* 
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Results  are  based  on  spectral  intensity  measured  at  7050  X.  Calculated  emission  included  continuum  only  for 
air  and  continuum  plus  Ha  line  for  Delrin.  See  Section  II,  part  D  for  discussion  of  I(x)  and  I*(x). 


SECTION  IX 


RESULTS  AND  CONCLUSIONS 

The  accomplishments  of  this  study  are  summarized  in  this  section.  Analytically  the 
Abel  inversion  has  been  extended  to  account  for  absorption  to  an  optical  depth  of  1  or  more 
and  applied  to  temperature  profile  determination;  information  has  been  obtained  on  the  pro¬ 
perties  of  Delrin  plasma  and  an  examination  has  been  made  of  the  significance  of  some  spectral 
lines  in  a  very  nonuniform  plasma. 

Experimentally,  a  stable,  high  pressure  plasma  was  generated  in  an  ablation  type  arc 
constrictor  and  the  arc  plasma,  under  conditions  of  the  experiments,  was  found  to  be  homo¬ 
geneous  along  the  arc  axis.  This  conclusion  is  based  on  the  luminous  homogeneity  of  the 
plasma  column,  negligible  pressure  rise  in  the  constrictor  during  arcing,  and  especially, 
the  constant  voltage  gradient  along  the  arc  axis.  The  voltage  ^gradient  could  be  measured 
much  more  precisely  than  in  previous  work  with  a  water  arc'’*'.  Extremely  high  gradients, 
(sometimes  exceeding  500  v/cm),  and  large  power  density  were  encountered.  In  view  of  its 
axial  homogeneity  the  ablation  type  constrictor  is  a  promising  tool  for  end-on  spectral 
measurements.  In  our  work,  side-on  measurements  required  a  quartz  window  which  caused 
some  difficulty  with  arc  stability,  though  not  enough  to  prevent  good  side-on  spectral 
measurements . 

A.  ARC  TEMPERATURE 

The  radial  distribution  of  the  arc  temperature  was  measured  for  250  ampere  arcs  in 
Delrin  plasma  with  pressures  in  the  100  to  150  atmosphere  range.  The  temperature  distri¬ 
butions  were  flat  with  an  occasional  slight  dip  at  the  axis  which,  in  view  of  the  difficulties  of 
both  the  analytical  and  experimental  methods,  was  assumed  to  be  spurious.  The  temperature 
measurements  were  based  on  the  theoretically  determined  intensities  at  the  wavelength  7050 
A.  These  measurements  depend  on  the  calculations  of  equilibrium  composition  (which  do  not 
include  the  possible  effect  of  lowering  of  the  ionization  potential  by  high  electron  density)  and 
absorption  coefficient,  as  discussed  in  Section  HI.  At  the  wavelength  chosen,  the  dominant 
continuum  radiation  is  precisely  that  which  is  best  known  (H^  and  bremsstrahlung).  The  Ha 
line  wing  is  also  well  known/*3) 

The  measured  temperatures  also  fall  within  a  range  bracketed  by  black  body  and  electrical 
conduction  considerations.  A  lower,  bound  can  be  deduced  from  the  black  body  limit  and  the 
fact  that  thermal  conduction,  as  deduced  from  the  measured  temperature  gradient  at  15,  000°K, 
and  the  thermal  conductivity  of  Yos^*),  is  less  than  10  percent  of  the  total  dissipation. 

Making  use  of  Figure  20,  the  radiation  from  the  current  carrying  part  of  the  arc  must  be 
more  than  90  kw/cm.  At  an  effective  diameter  of  0. 4  mm  radius,  this  corresponds  to  a  black  body 
temperature  of  16, 000°K  and  we  know  that  our  plasma,  while  not  optically  thin  is  not  black 
because  spectral  lines  can  still  be  seen  above  the  continuum  background.  On  the  other  hand, 
a  calculation  of  arc  current  using  the  measured  temperature  profiles  and  the  air  electrical 
conductivity  of  Viegas  and  Peng'23)  indicates  temperatures  are  not  excessive.  Calculated 
currents  are  already  about  50  percent  higher  than  measured,  both  in  Delrin  (Run  D-79)  and 
air  (Run  400). 
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B.  ENERGY  DISTRIBUTION 


The  distribution  of  continuum  radiation,  based  on  the  temperature  profile  of  Figure  31b 
as  computed  usingARCRAD  Ulb,  is  indicated  in  Figure  33.  The  plot  shows  the  relative  net 
radiant  energy  passing  radially  outward  at  any  radius.  Because  of  the  very  strong  continuum 
reabsorption  in  the  vacuum  ultraviolet,  60  percent  more  radiant  energy  passes  radius  0.4  mm 
than  eventually  emerges  from  the  arc,  and  most  of  what  does  emerge  is  in  the  visible  and 
near  UV.  This  again  demonstrates  that  the  vacuum  ultraviolet  is  a  powerful  mechanism  for 
dispersing  energy  in  the  arc  to  the  surrounding  gas. 

The  magnitude  of  the  calculated  radiation  leaving  the  current  carrying  part  of  the  arc  is 
30  kw/cm,  approximately  one-third  the  total  measured  dissipation.  This  is  low  but  may  be 
a  result  of  a  combination  of  line  radiation  and  some  unaccounted  for  continuum.  Most  of  the 
area  under  the  black  body  curve  at  these  temperatures  is  in  the  visible  and  near  UV  where 
the  radiation  processes  are  less  certain  than  in  the  infrared,  where  temperature  measure¬ 
ments  were  made. 

C.  SPECTRAL  LINES 

The  effect  of  lines  on  radiant  heat  flux  is  greater  than  was  anticipated  because,  rather 
than  representing  simply  a  narrow  band  of  black  body  radiation,  the  lines  are  so  greatly 
broadened  that  their  energy  is  distributed  over  a  wide  spectral  range.  The  Ha  line,  for 
example,  covers  about  1500  Two  of  the  lines  evaluated  using  the  ARCRAD  nib  program 
with  the  temperature  profile  of  Figure  31  are  shown  in  Figures  34  and  35.  Both  exhibit 
slight  self-reversal  which  could  not  be  confirmed  experimentally.  Energy  per  unit  arc 
length  for  the  Hq,  line  was  approximately  4  percent  of  total  calculated  continuum  and  for 
about  0. 6  percent. 

D.  APPLICATION  TO  AIR  ARCS 

The  motivation  for  this  work  was  a  need  for  information  about  arcs  in  air  at  high  pressure. 
However,  an  absolute  prerequisite  was  an  arc  sufficiently  stable  and  symmetric  to  yield 
quantitative  results.  Despite  intensive  efforts,  the  investigators  were  unable  to  stabilize  an 
arc  in  air  at  conditions  more  severe  than  about  100  amperes  at  100  atmospheres.  Therefore, 
use  was  made  of  ablation  stabilized  arcs  in  Delrin  plasma  with  the  results  being  applicable  to 
air  as  well.  Differences  in  transport  properties  are  small  and  the  radiation  is  approaching 
close  enough  to  black  body  to  preclude  large  differences  between  the  plasmas. 

Figure  32  shows  temperature  profiles  as  computed  from  continuum  radiation.  Aside 
from  uncertainties  already  discussed,  differences  are  a  result  of  the  combined  effects  of 
pressure,  constrictor  diameter,  current,  constrictor  type,  and  plasma  composition.  ■ 
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Pressure  has  a  relatively  small  effect,  as  can  be  seen  by  comparing  Runs  D-77  and 
D-79,  which  show  a  5  percent  change  in  centerline  temperature  for  a  50  percent  change  in 
pressure.  Runs  D-50  and  D-79  are  comparable  except  for  constrictor  diameter  (the  larger 
diameter  resulting  in  a  somewhat  broader  profile),  but  the  differences  are  small. 

There  is  a  greater  difference  between  the  250  ampere,  ablation  constrictor  Delrin 
plasma  arcs  and  the  air  arc*.  Though  no  direct  comparison  could  be  made  because  the  air 
arc  was  unstable  above  100  amperes  and  the  Delrin  arc  was  unstable  below  200  amperes,  it 
is  believed,  on  the  basis  of  indirect  evidence,  that  the  difference  in  the  profiles  is  primarily 
due  to  current  level.  Differences  in  transport  properties  are  small,  and  the  axial  uniformity 
of  the  ablation  type  arc  indicates  that  the  resulting  flow  field  has  little  effect  on  the  arc  core. 
Computed  continuum  emission  coefficients  at  X  =  7050  X  are  almost  identical  for  the  two 
plasmas  and  the  temperature. profile  is  not  sensitive  to  small  errors  in  i(r). 


*  Reevaluation  at  X  =  7050  A  has  resulted  in  a  lower  temperature  for  the  air  arc  than  previously 
reported^1). 
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CUMULATIVE  FRACTION  OF  CONTINUUM  RADIATION  LEAVING  ARC  COLUMN 


RADIUS  (MM) 

Figure  33.  Cumulative  Net  Continuum  Radiation,  Temperature  Profile  D-77 
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Figure  34.  Net  Radiant  Heat  Flux  Per  Unit  Arc  Column  Length  and  Per  Wave  Number 

in  Vicinity  of  O,^  Line 
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Figure  35.  Net  Radiant  Heat  Flux  Per  Unit  Arc  Column  Length  and  Per  Wave  Number  in 

Vicinity  of  Line 
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APPENDIX  I 


.ABEL  INVERSION  AND  INTEGRATION 
A.  ABEL  INVERSION  OF  A  POLYNOMIAL  FUNCTION 
Abel  inversion  is  the  evaluation  of  the  integral 

1  5  I*(x)dx 

f(r>  =  -i  /  (1) 

r 

and  is  the  inverse  of  Abel  integration  (Part  B).  Since  I(x)  is  a  measured  quantity  and  f(r) 
is  unknown,  the  inversion  process  is  the  key  one  in  the  temperature  profile  determination. 

It  depends  on  the  derivative  of  I(x)  so  that  the  manner  in  which  a  smooth  curve  is  fitted  to  the 
data  is  of  great  importance.  Cremer's  procedure  for  fitting  polynomials  by  the  method 
of  least  squares  has  been  used,  in  program  CF,  Appendix  n. 

The  arc  cylinder  is  divided  into  two  or  more  zones  and  a  separate  polynomial  is  computed 
for  each  zone.  For  each  polynomial,  all  the  data  within  the  zone  plus,  where  possible,  three 
data  points  beyond  both  inner  and  outer  zone  boundaries  are  used.  The  symmetry  require¬ 
ment  at  the  centerline  is  satisfied  by  using  only  even  order  terms  in  the  polynomial  associated 
with  the  innermost  zone.  The  polynomials  may  be  represented*  by 


I(x) 


m 

L 

3=1 


A.  x 
J 


2(j-D 


(innermost  zone) 


(2) 


I(x) 


m 

L 

3=1 


A.  x‘ 
3 


.3-1 


(other  zones) 


(3) 


*  Subscripting  is  compatible  with  computer  notation  which  forbids  zero  subscripts. 
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When  the  data  have  been  reduced  to  this  form,  the  Abel  inversion  may  be  carried  out 
analytically.  Differentiating  I(x)  and  integrating  term  by  term,  the  indefinite  integrals  are 
the  sum  of  a  final  series,  the  number  of  terms  depending  onthe  order  of  the  polynomials. 


Generalized  terms  could  be  written  for  Equations  4  and  5  but  they  are  cumbersome 
and  the  pattern  for  extension  to  higher  order  is  clear. 
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Note  that  the  integral  is  evaluated  at  a  fixed  value  of  r.  which  is  also  the  first  lower 
limit  of  x.  Successive  upper  and  lower  limit  values  of  x  occur  at  the  zone  boundaries, 
where  the  array  of  coefficients  changes,  reaching  a  final  upper  limit  at  the  arc  radius. 
Evaluation  of  Equations  4  and  5  is  performed  by  subroutine  THIN  of  the  optically  thick 
inversion  program,  OTHICK,  Appendix  II. 

B.  ABEL  INTEGRATION  OF  A  FUNCTION  KNOWN  AT  DISCRETE  POINTS 

Two  different  Abel  integrations  are  required  in  the  non-optically  thick  inversion  process, 
one  to  find  the  intensity  profile  I*(x),  .and  one  to  determine  the  optical  depth  of  the  absorbing 
path  as  a  function  of  physical  location.  The  Abel  integral  can  be  thought  of  as  a  summation 
over  the  volume  elements  along  a  line  of  sight,  off  axis  in  general,  through  the  arc  column 
Figure  36.  If  absorption  is  not  significant  the  integral  is  given  by 

V  -  ■/» 


The  factor  of  2  arises  because  the  line  of  sight  is  twice  the  limits  of  integration.  When  absorp¬ 
tion  is  significant,  emission  from  the  volume  elements  behind  the  centerline  will  be  asorbed 
to  a  greater  extent  than  from  those  in  front  of  it  and  instead  of  the  factor  of  2  there  will  be  two 
separate  integrals,  with  differing  f(r). 

In  our  situation,  the  function  f{r)  has  been  determined  by  Abel  inversion  at  uniformly 

spaced,  discrete  points,  r^  . . .  r^,  r.+^. . .  r  ,  and  x,  the  lower  limit  of  integration,  is  equal 

to  r^,  see  Figure  36.  Integration over  each  interval  r^  to  r.  is  performed  with  several 

variations  on  the  theme  of  Simpson's  rule,  i.  e. ,  by  approximating  the  true  curve  with  a 

parabola  through  f.  ,,  f.,  and  f.  When  r.  ,  >r,  we  let 
i-1  l  x+1  i-1  k 


and,  with  the  notation  f  -  f  (rj),  find  that 


61.  = 
i 


(ri+l~ri^ 

12 


(-f.  +  8f.  +  5f  ) 

'  l-l  l  i+l' 


(k  >  1;  i  s  k  +  2) 


(8) 


When  either  r ^  =  r  or  r-_^  =  r  Equation  7  blows  up.  The  integral  remains  finite  however, 
so  Equation  6  was  rewritten,  defining  rj  =  r/r^  as 
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I(*7)  a  r. 


V 

f  MT7 )  V  d  r? 

•/  Vr? 2  -  1 


(9) 


Integration  by  parts  yields , 


im 


I  =  rk  [f  (T7)  V772  -  1  -  /  f»  (77)  Vt?2  -  1  dT?]1 


(10) 


Performing  the  indicated  operations  with  the  parabola  through  1  f,  f  results  in 


SIi =  \  {fm  Vv/ - 1  - fi 


-  ir  fvi  Vvi4  - 1  -  ^ 


1  -  An 


vS:1 

,”1 +V”i2  - 1 


0, 


r  3/2 

n  2  1)  ( 

1  ,  ] 

lyi+i  V  1 

>  Ji 

=  k  or  i  =  k  +  1)  (11) 


7?.  =  r./r 
x  1  k 


(12) 


0-.  = 


[f..i  -f.  .3 


1  2  <Vi  “V  i+1  i-1 


(13) 


0  =  ■ 
2 


Cf.  ,  "2f  +  f  ] 


(Vi"V 


2  "i-l  i  i+1 


(14) 
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(15) 


with  6  I.  given  by  Equations  8,  11,  16,  or  17  as  appropriate.  The  integration  process  is 
performed  by  subroutine  INT3P  in  the  optically  thick  inversion  program,  INVER. 
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Figure  36.  Abel  Integration 


APPENDIX  II 


NON-OPTICALLY  THIN  ABEL  INVERSION  PROGRAMS 

A.  INTRODUCTION 

The  computer  programs  CF  and  INVER,  developed  to  convert  raw  intensity  data  to  tem¬ 
perature  profiles,  are  discussed  in  this  Appendix.  Table  IV  presents  a  brief  summary  of 
their  purpose,  inputs  and  outputs.  Computations  were  carried  out  on  the  GE  Desk  Side  Com¬ 
puter  System  (DSCS  n),  which  is  a  time  sharing  system  using  a  standard  teletypewriter  as  a 
remote  terminal.  The  system  uses  FORTRAN  IV  with  a  number  of  special  features.  Those 
special  features  and  those  library  subroutines  and  built-in  functions  which  pertain  to  CF  and 
INVER  are  summarized  in  Tables  V  and  VI. 

One  of  the  special  features  is  a  disc  storage  capability.  Outputs  of  program  CF  are 
stored  automatically  in  binary  disc  files  and  are  retrieved  as  desired  when  program  INVER 
is  run.  This  Is  a  convenience  but  not  a  necessity,  and  other  means  of  data  transfer  can  be 
employed  with  minor  modifications  to  the  programs. 

B.  PROGRAM  CF 

Program  CF  fits  polynominals  to  raw  intensity  measurements  prior  to  Abel  inversion. 
Data  may  be  at  random  intervals  and  both  the  degree  of  the  polynominals  and  the  number  of 
zones  (i.e. ,  number  of  different  polynominals)  can  be  varied  to  obtain  best  fit.  This  a  mat¬ 
ter  of  some  judgment,  particularly  since  it  is  the  slope  of  the  intensity  curve  that  is  of  inter¬ 
est,  but,  as  a  guide,  program  CF  computes  root  mean  square  deviation,  cr,  for  each  fitted 
curve. 

If  the  data  are  "reasonable,  "  o  wiU  tend  to  decrease  as  degree  is  increased*  until  a 
minimum  is  reached.  Further  increases  in  degree  then  cause  c  to  remain  roughly  constant, 
sometimes  increasing  or  oscillating  slightly.  When  this  happens  the  additional  terms  in  the 
polynominal  are  simply  putting  extra  "wiggles"  in  the  curve  and  the  result  is  a  progressively 
poorer  representation  of  the  "true"  slope.  Third  degree  consistently  gave  best  results  with 
our  data,  though  the  program  arrays  are  sized  for  up  to  fifth  degree.  At  least  two  zones 
must  be  used  because  of  the  special  treatment  given  the  center  zone  to  assure  zero  slope  at 
the  centerline  (Appendix  I).  The  program  variables  are  sized  for  up  to  five  zones,  and  a 
total  of  50  data  points. 

Data  are  read-in  in  dimensional  form  with  the  units  specified  in  Table  IV,  but  are  normal¬ 
ized  prior  to  computation  of  the  least  squares  fit  to  avoid  exponential  overflow  or  underflow. 
Also,  each  curve  is  fitted  to  its  own  origin  of  coordinates  to  minimize  round  off  error  in  the 
least  squares  algorithm.  Polynominal  coefficients  are  computed  using  library  subroutine 


*  Polynominals  must  be  at  least  third  degree  for  proper  operation  of  INVER. 
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LINEQ,  according  to  the  least  squares  criterion,  as  discussed  by  Scarboro^).  The  normal¬ 
ization  process  is  then  reversed  and  the  coefficients  transformed  to  a  common  axis.  The 
latter  process  requires  that  a  polynomial  of  the  form* 

y  =  \  +  ^2  (x  -  a)  +  Ag  (x  -  a)2  +  .  .  .  .  +  An+1  (x  -  a)n  (1)  . 

be  transformed  to 

y  =  A1  +  A2x  +  A3x2  +  .  ..  .  +An+iXn  (2)  . 

The  process  involves  finding  the  A's  as  functions  of  (a,  A . by  expanding 

Equation  1  and  equating  like  powers  of  x.  Library  subroutine  BICOF  was  used  for  generation 
of  the  necessary  arrays  of  binomial  coefficients. 

C.  PROGRAM  INVER 

Program  INVER  is  used  to  compute  temperature  profiles.  It  has  as  its  input  the  poly¬ 
nomials  representing  observed  intensity,  another  polynomial  relating  emission  coefficient 
to  temperature,  and  certain  run-identification  information.  Because  emission  coefficient 
i(r)  can  best  be  represented  as  a  polynomial  in  temperature  rather  than  vice  versa,  Newton's 
method  for  finding  the  root  of  a  polynomial  was  used  to  obtain  T  for  a  given  i(r).  Thus  it  is 
vital  that.the  emission  coefficient  polynomial  increase  monotonically  with  temperature  in  the" 
range  of  interest.  A  curve  fit  procedure  was  used  which  evaluated  all  roots,  extrema,  and 
points  of  inflection  so  that  complete  behavior  could  be  ascertained  (see  sketch). 
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Choice  of  the  number  of  discrete  points  for  numerical  integration  is  governed  by  a  bal¬ 
ance  among  truncation  and  round-off  errors  and  computation  time,  but  since  all  the  numerical 
processes  are  integrations  (differentiation  of  intensity  profiles  is  analytic.  Appendix  I), 
computational  errors  are  not  a  serious  problem,  and  15  to  20  points  on  the  profile  have 
proved  adequate. 

The  block  diagram  of  Figure  37  shows  the  steps  in  the  computation.  From  one  to  six 
.iterations  were  necessary  to  converge  to  minimum  rms  deviation,  between  measured 

and  computed  I(x).  Typically,  am£n  was  between  1  and  2  percent  of  1(0). 

The  following  subprograms' form  an  important  part  of  INVER: 

Function  ROOT  is  a  separate  function  subprogram  for  calculating  -  v2,  which 
is  useful  in  evaluating  both  the  Abel  integral  and  its  inversion. 

Subroutine  IOFX  is  a  straightforward  calculation  of  I*(x)  at  discrete  points  from 
the  polynomials. 

Subroutine  THIN  performs  the  Abel  inversion  of  I*(x).  The  inversion  follows  the 
analytic  procedure  outlined  in  Appendix  I  and  is  performed  at  each  successive 
Xfc  =  r^  from  centerline  to  outside  edge.  The  indicator  LIM  is  set  as  follows: 

LIM  Significance 

f 

1  Lower  limit  of  integration  for  zones  outside  the  zone  for  which 
LIM  =  3 

2  Upper  limit  of  integration  for  zones  outside  the  zone  for  which 
LIM  =  3 

3  Upper  and  lower  limits  of  integration  at  the  first  zone  boundary 
where  the  radius  exceeds  r^ 

LIM  =  3  is  a  special  case  because  then  the  lower  limit  is  zero  except  for  the  logarithm  term 
.  in  Equation  5  of  Appendix  I,  and  this  is  included  by  setting  DEN  =*  x^  *  r^. 

Subroutine  POLYTEM  is  a  straightforward  evaluation  of  i(r)  and  its  first  derivative 
(needed  for  Newton's  method  correction  of  trial  value  of  temperature)  at  a  specified 
temperature. 

Subroutine  INT3P  performs  the  Abel  integration  of  f(r)  following  the  procedure 
described  in  Appendix  I. 

<# 

Subroutine  FITX  fits  new  polynomials  to  corrected  values  of  I*(x)  in  the  same  way 
.  as  the  original  data  prior  to  Abel  inversion.  Now,  however,  the  zone  structure  and 

order  are  already  established  and  the  points  are  at  uniform  intervals  and  can  be 
expected  to  conform  very  closely  to  a  smooth  curve. 
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Table  IV.  Programs  for  Computation  of  Temperature  Profiles 


Program 

Purpose 

Input 

Output 

CF 

Converts  raw  intensity  data  to 

1. 

Material  (e.g.  air,  Delrin)  pressure, 

1. 

(optional)  Column  listing  of  input  2 

polynomials  and  transmits 

wave  number  of  intensity  measure¬ 

other  Input  information  to 

ments  and  coefficients  of  the  poly¬ 

2. 

Root  mean  square  deviation  for  each  zone 

program  INVER. 

nomial  representing0cmission  coeffi¬ 
cient  (erg/cm^  sec-A-str) 

3. 

(optional)  Table  of  values  for  zone: 
displacement,  raw  data,  smoothed  data, 

2. 

Number  of  data  points,  displacement 
from  axis  (cm)  of  each  data  point. 

differences 

intensity  measurements  (10^  erg/sec- 
cm  -X-str).  Data  must  be  in 

4. 

Items  1  and  3  of  input 

order  of  increasing  displacement 
from  the  axis 

5. 

Coefficients  of  the  polynomials 

3. 

Degree  of  polynomials  (3  to  5), 
number  of  zones  (2  to  5),  arc 
radius  (cm) 

INVER 

Performs  Abel  inversion 

1. 

Number  of  points  to  be  computed 

1. 

Average  difference  and  root  mean 

accounting  for  effects  of  self 

for  each  zone 

square  deviation  between  measured  I(x) 

absorption 

2. 

Output  from  program  CF 

and  computed  I(x)  for  each  iteration 
(see  Section  IID). 

3. 

Integer  Index  after  each  iteration 

2. 

Table  of  intensity,  emission  coefficients. 

(0  =  repeat  cycle,  1  -  print  output, 

and  temperature  as  functions  of  radius 

then  repeat  cycle,  2  =  new  case) 

(see  input  3) 

Table  V.  Special  Features  of  the  FORTRAN  Language  Used  for  Programs  CF  and  INVER 
Symbols : 

*  -  Comment 

&  -  Continuation  of  statement 

Statements; 

ASC  n 

DATA 

FILENAME 

PRINT: 

READ; 

Other  Notes: 

Mixed  mode  expressions  are  allowed  and  were  used. 

Logical  IF  is  used  in  CF . 

Subscripts  of  an  array  need  not  be  shown  explicitly  in  arguments  of  subroutine. 


Specifies  alphameric  variable  (used  to  record  the  type  of  plasma, 
e.  g. ,  Delrin  or  air) 

Statement  enabling  one  to  put  many  program  constants  into  a 
single  statement 

Used  to  specify  the  name  of  the  disc  storage  file  (used  for  transfer 
of  CF  output  to  INVER).  STOW  is  a  filename  variable,  whose 
current  "value"  is  the  name  assigned  to  the  binary  disc  file  con¬ 
taining  data  from  a  specific  run.  The  experimental  run  designation 
was  used  to  name  the  file. 

Permits  printout  of  alphmeric  information  without  need  of  a 
FORMAT  statement 

Used  for  free  format  data  input 
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Table'  VI.  ‘Library  Subroutines  and  Built-in  Functions  Used  by  Programs  CF  and  INVER 
LINEQ  (AM,  AB,  NA,  NB,  IDIM) 


Solves  an  array  of  simultaneous  linear  equations  by  Gaussian  elimination 

AM  Array  containing  the  matrix  of  coefficients 

AB  Array  containing  right-hand-side  vector  on  input  and 

(  solution  vector  when  completed 


NA 


f,r  i 

NB  ' 


IDIM 

.  ^  * 


Number  of  equations 

Number  of  right -hand-side  vectors 

First  dimension  of  AM  and  AB  arrays 


BICOF  (0/  N,‘  COF) 

Generates  an  array  of  binomial  coefficients 

COF  (J,l)  Array  containing  the  coefficients  (J  =  1,  N) 

N  Degree  of  the  array.  For  example  if  N=4  the  array  would 

be  (1,4, 6, 4,1) 


Built  in  Functions  Used  by  CF  and  INVER 


EXP 

Exponential  function 

ALOG 

.  /  1 

Natural  logarithm 

SQRT 

Square  root 

ABS 

Absolute  value 
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D.  GLOSSARY 


AB  (J,  1) 

ABS  (K) 

t  <  * .  r  ■ 

.  •  * 

AJ  (J) 

AM  (JK,  JL) 
BICOF 
C  (J,  1)  j 

c  (J,  2)’  [  .  ... 

CF 

COF  (JR,  I).;' 

COF  (JR  J+l) 

■  / 1  - 

CLG 

CP 

CW5  , 

C2W 

'  ,  »  »  %'•  i  *  . 

DE 

DEN 

DI 

DIS 

DIT 

DR 

DT  (K) 


Array  for  LINEQ  (see  Table  VI) 

Local  absorption  coefficient 

Coefficients  of  polynomial  for  emission  coefficient  as  f  (T) 

Array  of  coefficients  of  the  polynomials  for  I*(x) 

See  Table  IV 

Constants  for  calculation  of  coefficients  of  terms  in  Abel  inversion,  see 
Equations  4  and  5,  Appendix  I 

Curve  fit  program 

Output  of  BICOF  (Table  VI) 

Modified  binomial  coefficient  array  for  translation  to  common  axis  (see 
part  B  of  this  Appendix) 

Coefficient  of  logarithmic  term  in  Abel  inversion,  Equation  5,  Appendix  I 
Coefficient  of  each  term  in  Abel  inversion,  Equations  4  and  5,  Appendix  I 
Planck  function  term 
Exponent  in  Planck  function 
Increment  in  E 

Lower  integration  limit  value  of  argument  of  logarithm  term  in  Abel 
inversion,  Equation  5,  Appendix  I.  See  Subroutine  THIN,  part  C  of  this 
Appendix. 

Difference  between  calculated  and  observed  I(x) 

Running  sum  of  DI  squared;  eventually,  rms  deviation 
Running  sum  of  DI;  eventually,  mean  deviation 
Radius  increment 

Incremental  Abel  integral,  Subroutine  INT3P 
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DTIS  (KX) 


DT  applied  to  I*(x) 


DTXR  (KX) 

E 

ER 

FITX 

INT3P 

INVER 

10  FX 

J 

JK  ) 

JL  j 

JM 

JORD 

JP 

K 

KDIF 

KF 

KFIN 

KL 

KM 


DT  applied  to  optical  depth  for  computation  of  Cl(x)]cajc;  also  used  by 
Subroutine  POLYTEM  for  derivative  of  temperature-emission  coefficient 
function 

Eta,  dimensionless  argument  in  Abel  integral,  Equation  10,  Appendix  I 

Vl  - 1 )2,  Equation  10  Appendix  I 

Subroutine,  see  part  C  of  this  Appendix 

Subroutine,  see  part  C  of  this  Appendix 

Non-optically  thin  Abel  inversion  program 

Subroutine,  see  part  C  of  this  Appendix 

Index  of  intensity  polynomial  coefficients  and  emission  polynomial 
coefficients 

Indices  in  least  squares  algorithm 
Maximum  value  of  J 

Degree  of  emission  coefficient  polynomial 

J-l,  argument  of  BICOF 

Index  on  radius  or  displacement  from  axis 

KSTRT-KI 

Largest  value  of  K  for  which  RD  (K)  falls  at  or  within  the  boundary  of  a 
particular  zone 

Largest  value  of  K  for  set  of  data  to  which  a  polynomial  will  be  fitted, 
see  KSTRT 

Smallest  value  of  K  for  which  RD  (K)  falls  at  or  within  the  boundary  of  a 
particular  zone 

Maximum  value  of  K 
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KMX 


KM-1 


KSTRT,.  ? 


KX 

KX1 

KX2 

KZ 

L 

LIM 

UNEQ 

M 

MAT 

MDATA 

MS 

MZ  (NZ) 
N 

NA,  NB 
NNEW 

i  -  •  i 

NNN 

NNZ 

NORD 

NP 


First  value  of  K  for  which  RD  (K)  falls  within  the  set  of  data  to  which  a 
polynomial  will  be  fitted  (up  to  3  points  beyond  zone  boundary  are 
included  where  possible) 

Summation  index;  KX  ^  K 

KX  +1 

Index  on  first  point  of  3 -point  Abel  integral,  KX2  =  2  or  KX+2 

Number  of  subintervals  per  zone  (INVER  only). 

Counts  number  of  iterations.  Initial  (optically  thin)  inversion  is 
numbered  zero. 

See  Subroutine  OTHIN,  part  C  of  this  Appendix. 

See  Table  VI. 

Number  of  data  points  included  in  a  curve  fit  calculation 

ASC  II  variable;  Table  VI 

Total  number  of  data  points 

Indexing  constant,  Subroutine  IFTX 

Value  of  K  at  outer  boundaries  of  zones 

NORD  +  1 

See  UNEQ 

Option  indicator.  Used  by  operator  to  execute  choices  to  be  made  while 
program  is  being  run. 

See  NNEW 

Total  number  of  zones  (2  s  NNZ  ^  5) 

Degree  of  I(x)  polynomials 
See  NNEW 
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NPO 

NXY 

NZ 

P 

PHI  (I,  J) 

PHI) 

PH2 
PH3  ) 

PI 

POLYTEM 

PR 

PTSF 

Q  (KX) 

QI(KX) 

R(K) 

RD  (K) 
RDMAX  (NZ) 
RDMIN  (NZ) 
REF 
ROOT 
RO 


Counts  iterations  to  find  i(r)  from  T(r).  Causes  print  out  if  this  process 
does  not  converge. 

See  NNEW 

Index  on  zone 

Pressure  (program  CF) 

Array  of  values  of  polynomial  functions  for  least  squares  alogorithm. 
Coefficients  of  Abel  integral.  Equation  12,  Appendix  I 

3.1415927 

Subroutine,  see  part  C  of  this  Appendix 
Pressure  (program  INVER) 

Abel  integral  function  defined  at  beginning  of  subroutine  INT3P;  see 
Equation  11  of  Appendix  I 

Array  of  inputs  to  Abel  integral,  dummy  argument  of  the  subroutine 
INT3P 

Same  as  Q  to  start  with,  but  modified  by  Subroutine  INT3P.  Used  to 
avoid  changes  in  the  array  in  INVER  which  Q  (kx)  represents 

Radius  (program  INVER) 

Radius  (program  CF) 

Radius  of  outer  boundary  of  zone 

Radius  of  inner  boundary  of  zone 

Reference  value  for  normalization  prior  to  curve  fit 

Function,  see  part  C  of  this  Appendix 

Arc  radius 
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SIGMA 

STOW 
T  (K) 
THIN 

lO'Ji 

TI  (K) 

TIA  (K) 
TICK 
TILG 
TIS  (K) 

TIP 

TLX  (K) 

.  v> 

TIZ 

TPW  ■ 
TXR 

WN 

X 

X(I) 

YA  (K) 
YI(K) 


rms  deviation:  raw  data  from  polynomials  in  program  CF;  computed 
intensity  from  measured  intensity  in  program  INVER 

See  Table  V 

Temperature 

Subroutine ,  see  part  C  of  this  Appendix 

Emission  coefficient 

Intensity  including  self-absorption 

Residual  in  calculation  of  temperature  from  emission  coefficient 

Coefficient  of  logarithm  term  in  Abel  inversion,  Equation  5,  Appendix  I 

I*(x);  intensity  which  would  have  been  observed  if  part  of  the  energy  had 
not  been  reabsorbed. 

Running  sum  in  Abel  inversion 

Observed  intensities,  as  computed  from  polynomials  fitted  to  data 
Running  sum  to  compute  TI  (K) 

Running  sum  to  compute  T  (K) 

Optical  depth;  also  used  by  Subroutine  POLYTEM  for  temperature 
calculation 

Wave  number 

Displacement  from  arc  centerline;  this  the  lower  limit  of  Abel  integral 
(program  INVER) 

Dependent  variable,  least  squares  algorithm 
Intensity  from  fitted  polynomials  (program  CF) 

Intensity  data  (program  CF) 
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E.  LISTINGS 


00010  * 
00020 
00030 
00040 
00050 
00060 
00070 
00080 
00090 
00100 
00110 
00120 
00130  * 
00140  * 
00150 
00160 
00170 
00180 
00190 
00200 
00210 
00220 
00230 
00240 
00250 
00260 
00270 
00280 
00290 
00300 
00310  * 
00320  « 
00330 
00340 
00350 
00360 
00370 
00380 
00390 
00400 
00410 
00420 
00430 
00440 
00450 
00460 
00470  * 
00480  * 
00490 
00500 
00510 
00520 
00530 
00540 


CF#  UP  TO  50  DATA  POINTS*  5  ZONES*  5TH  DEGREE  POLYNOMIALS 
1  FORMATC  1H  •  F6.  3*  1XE1 1 . 3> 

3  FORMATC  "  "> 

8  FORMATC  3C 1PE1 1 • 3) *  2C 1PE10*2) ) 

10  F0RMATC2C 1PE10.2)) 

COMMON  ABC  6*  1  )  *  ARC  6*  6)  •  COFC  6*  7> 

COMMON  XC20)  *YC20)*PHIC50*6)*AMC  6*  6> 

COMMON  RDMAXC  6) *  RDMINC  6) *  RDC  50) *  YI C 50) 

COMMON  YAC16) 

DIMENSION  AJC  7) 

FILENAME  STOW 
ASCII  MATC1) 

READ  IN  DATA 
100  PRINT  3 
PRINT  3 

PRINT* " AI R  OR  DELR*  PRESS*  WN»  JORD*  JM*  AJCJ)" 

PRINT  3 

READ*  MAT*  P*  WN*  JORD*  JM*  CAJCJ)*  J=1.JM) 

PRINT  3 

PRINT*  "MDATA*  ALL  RCCM).  ALL  I C X ) C 10. E6ERG/SEC-CM2- A- STR) " 
PRINT  3 

READ *M DATA* C  RDC  K ) »K“1 *MDATA) *CYICK)»K»1 #  MDATA) 

RDCMDATA+ 1 )» RDC M DATA) ♦ 1 • 

DO  102  K»l* MDATA 
102  YICK)=YI CK)*1.E6 

PRINT* "LI  ST  DATA*  1  *  " 

READtNNN 

IFCNNN)  105*105*104 

104  PRINT  1  *  C  RDCK) » YI CK) »K* 1* MDATA) 

SET  UP  ZONE  STRUCTURE 

105  PRINT* "NO RD*NNZ * R0" 

READ*  NORD*  NNZ*  R0 
N«NORD+ 1 

RDMINC l)a0. 

RDMAXC  1)*R0/NNZ 
DO  110  NZ*2*  NNZ 
RDMINCNZ)*RDMAXCNZ- 1 ) 

110  RDMAXCNZ)*RDMAXC 1>*NZ 
RDMAXCNNZ  )*R0 
PRINT*”PRINT  COORDs  1*  •• 

READ*  NXY 

IFCNXY)  1000*  114*  112 
112  PRINT*"  R  YI  YA” 

114  DO  300  NZ* 1  *  NNZ 

SELECT  DATA 
K*0 

120  K*K+ 1 

I FC  RDCK ) *  RDMAXC NZ ) - .0000 1)1 20* 1 20*  130 
130  KFIN*K 
K*K+2 

I FCMDATA-K)  150*  140*  140 
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00550 

140 

M*K 

00560 

GO  TO  160 

00570 

150 

M«MDATA 

00580 

KFINsMDATA 

00  590 

160 

K*0 

00600 

170 

K*K+  1 

00610 

I FC  RDCK  > -RDMIN  CNZ )♦• 0000 1 )  170* 180# 180 

00620 

180 

KSTRT=K - 1 

00630 

KI=K-3 

00640 

IFCKI)  190#  190#  200 

00650 

190 

KSTRT*  1 

00660 

KI*  1 

00670 

200 

KDIF*KSTRT-KI 

00680 

M*M+1-KI 

00690 

DO  210  1*1#  M 

00700 

XCI >*RDCI-1+KI ) 

00710 

210 

YCI )  *YI Cl- 1*KI ) 

00720 

* 

00730 

♦ 

LEAST  SQUARES 

00740 

REF*XCM)-XC 1 > 

00750 

DO  255  I*  1  #M 

00760 

PHICI# 1)*1. 

00770 

IFCNZ-1)  1000#  220#  240 

00780 

220 

DO  230  J*  2#  N 

00  790 

230 

PHI  Cl# J)*CXCI >/REF)**C2*CJ- 1>> 

008  00 

GO  TO  255 

00810 

240 

DO  250  J*2#  N 

00820 

250 

PHI CI,J)*CCXCI)-XC 1)>/REF)**CJ-1) 

00830 

255 

CONTINUE 

00840 

DO  256  JL«1#N 

008  50 

DO  256  JKa 1 #N 

008  60 

AMCJK# JL)«0. 

00870 

DO  256  I-1«M 

00880 

256 

AMCJK# JL)=AM<JK#JL)+PHI (I# JL>*PHI (I# JK) 

00890 

DO  257  JK* 1 #N 

00900 

ABC JK# 1 )*0. 

00910 

DO  257  I* 1 #M 

00920 

257 

ABC JK# 1 )■ ABC JK# 1) +YCI>* PHICI# JK) 

009  30 

CALL  LINEQC AM#  AB#N# l#  6> 

009  40 

DO  261  J*2#N 

009  50 

IFCNZ-1)  1000#  258#  260 

009  60 

258 

ABC J# 1 )*ABC J# 1 ) * C  REF**C  2* C 1-J>>> 

009  70 

GO  TO  261 

00980 

260 

ABC  J# 1 ) * ABC  J# 1)*C REF**  C 1-J) ) 

0099  0 

261 

CONTINUE 

01000 

* 

01010 

* 

SHIFT  AXIS 

01020 

1FCCNZ.EQ.  D.OR.  CXC  1>.EQ.0.  >>  GO  TO  266 

01030 

DO  262  J* 1 # N 

01040 

JP=J-1 

01050 

CALL  BICOFC0# JP#COF> 

01060 

DO  262  JR*  1#  J 

01070 

262 

COFC  JR#  J*  1  )*COFC  JR#  1  > 

01080 

DO  264  JR* 1 #  N 
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01090 
01100 
01110 
01120 
01130 
01140 
01150 
01160 
01170 
01180 
01190 
01200 
01210 
01220 
01230 
01240 
01250 
01260 
01270 
01280 
01290 
01300 
01310 
01320 
01330 
01340 
01350 
01360 
01370 
01380 
01390 
01400 
01410 
01420 
01430 
01440 
01450 
01460 
01470 
01480 
01490 
01500 
01510 
01520 
01530 
01 540 
01550 
01560 
01570 
01580 
01590 
01600 
01610 
01620 


ARC  JR*NZ)*0« 

00  264  J*  JR*  N 

264  ARC  JR*NZ)b  ARC  JR*NZ)  +  C  C-1>**CJ-JR> >  * ABC J* 1>* 

6  CRDCKI>**CJ-JR))*C0FCJR*J+1> 

GO  TO  270 
266  DO  268  J«1*N 
268  ARC J*NZ )* ABC J*  1) 

* 

*  CALC  AND  PRINT  COORD 
270  M»KDI  F+KFIN-KSTRT+ 1 

I  STRT*KDI  F+ 1 
SIGMA«0. 

DO  290  IalSTRT*  M 
IFCNZ-l)  1000*  272*  276 
272  CONTINUE 

YAC I  )■  ARC  1*1) 

DO  274  J*2*N 

274  YACI)aYACI >+ARCJ* 1  )«XC I )«*C 2*C J- 1 > > 

GO  TO  290 
276  CONTINUE 

YAC  I  )- ARC  1*NZ)+ARC2*NZ)*X(  I  ) 

DO  280  J*3*N 

280  YAC I >»YACI)+ARCJ.NZ)*CXCI >**CJ-1>> 

290  SI  GMA«SI  GMA+C  YC  I  >  -  YAC  I  >  >**2 
SI  GMA» SORT  CSIGMA/CM-KDIF)  > 

IFCNXY)  298*298*294 

294  CONTINUE 

DO  29  6  1*1  STRT*M 

295  PRINT  8*  XC I ) * YC I > * YAC I ) 

29  6  CONTINUE 

298  PRINT  10*  SIGMA 
PRINT  3 
300  CONTINUE 

* 

*  RESULTS  TO  FILE  OR  TRY  AGAIN 

PRINT*  "NEW  FI  T**0*  TO  FILE«1*  " 

READ*  NNEW 

IF  CNNEW)  1000*  105*  302 
302  PRINT* "FILENAME" 

READS  STOW 
BEGIN  FILE  STOW 

WRI  TEC  STOW)  MAT*P*  WN»  JORD*  JM*  C  AJC  J)  *  J*  1*  JM)  *N*NNZ*  R0 
END  FILE  STOW 
DO  310  NZ* 1 *  NNZ 
WRI  TECST0W)CARCJ*NZ),J«1*N> 

310  PRINT*  CARCJ*NZ>*  J«1*N> 

PRINT  3 
END  FILE  STOW 

PRINT*  "NEW  FIT=0*  NEW  DAT  As  1" 

1  CLOSE  FILE  STOW 
READS  NNEW 

IFCNNEW)  1000*  105*  100 
1000  STOP 
END 
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00010 

00020 

00030 

00040 

000S0 

00060 

00070 

00080 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

00160 

00170 

00180 

00190 

00200 

00210 

00220 

00230 

00240 

00250 

00260 

00270 

00280 

00290 

00300 

00310 

00320 

00330 

00340 

00350 

00360 

00370 

00380 

00390 

00400 

00410 

00420 

00430 

00440 

00450 

00460 

00470 

00480 

00490 

00500 


■  INVER*  2ND  TO  5TH  ORDER*  2  TO  5  ZONES*  UP  TO  10  PTS/ZONE 

COMMON  ARC  6*  5>  *  MZC  5) *  CC6*2>*  R<5!>*  TIS<51>»  TIXC  51 > 
COMMON  ABSC51 ) *  TC 51 ) *  TIC5D*  TIAC51 > 

COMMON  TXRCSI )*  DTISC5D*  DTXRC51 )*  DR*X*K*KX1 *KMX*KM 
COMMON  PHIC17*6)*AM(6*6>*C0FC6*7)*YC17)*ABC6*1 ) 

COMMON  AJ<7>*  JORD* JM 
FILENAME  STOW*  P 
ASCII  MATCH 

1  FORMAT.C  4C1  PE  11*31) 

2  FORMATC 1H  *A4* F6. 0* F8.0* 13*15* 14) 

3  FORMATC"  "> 

i 

►  INPUT*  ZONE  STRUCTURE  AND  PROGRAM  CONSTANTS 

PI »3. 141 5927 

100  PRINT* "FILENAME*  PTS/ZONE  " 

READ*  STOW*  KZ 
BEGIN  FILE  STOW 

READC  STOW* END= 1 0S >MAT* PR*  WN* JORD* JM*  C A JC J> *  J= 1 • JM) * N*NNZ* R0 
105  DO  110  NZ= 1 *NNZ 

READC  STOW*  END= 110) C ARC J*  NZ  > • J» 1 • N) 

110  MZCNZ)«KZ*NZ+1 
.  KM«MZCNNZ) 

KMX*KM-1 
DR-R0/KMX 
DO  120  K»1  *  KM 
120  RCK)»CK-1)*DR 

►  CW5*CERG/SEC-CM2-A-STR) 

CW5«CWN**5)*1.1909E-13 
C2W=1 .438*WN 
C<  2* 1 )*2» 

DO  130  J&3*  N 

CCJ*1 )=C2.*C J-i ) )/C2.*J-3« ) 

130  CC J*2)»C J-l*)/CJ-2.) 

¥ 

¥  ISTARCX)  TO  ICR)  AND  TEMPERATURE 
TI CKM )«0* 

TlACKM)*0* 

CALL  IOFXCKZ*  N*NNZ) 

&  y ■;  DO  150  K*l*  KM 
150  TIXCK)t*TISCK) 

NORDsN- 1 
PRINT  3 

PRINT*"  PC ATM)  WN  JORD  NORD  NNZ" 

PRINT  2*  MAT*  PR*  WN*  JORD*  NORD*  NNZ 
L«-l 

220  L*L+1 

CALL  THINCN*NNZ) 

DO  570  Ke 1 «KMX 

T!CK)a-TICK)/PI 

TCK)al0000. 
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00510 
00520 
00530 
00540 
00550 
00S60 
00570 
00580' 
00590 
00600 
00610 
00620 
00630 
00640 
00650 
00660 
00670 
00680 
00690  * 
00700  * 
00710 
00720 
00730 
00740 
00750 
00760 
00770 
00780 
00790 
00800 
00810  * 
00820  * 
00830  * 
00840  * 
00850 
00860 
00870 
00880 
00890 
00900 
00910 
00920 
00930 
00940 
00950 
00960 
00970  * 
00980  ♦ 
00990 
01000 


CALL  POLYTEM 

IFCTI CK)-TXRCK> )  570#  570#  300 
300  TCK)**TCK)  42000# 

IFCTCK)- 30000)  360#  360#  340 
340  PRINTfTEMP  TOO  HIGH" 

360  NPO»0 

CALL  POLYTEM 

IFCTI CK)-TXRCK) )  400#  570#  300 
400  TCK)«TCK)-CTXRCK)-TICK))/DTXRCK) 
NPOaNPO+1 

IFCNPO-50)  430#  430#  420 
420  PRINTs"T  UNSTABLE" 

430  CALL  POLYTEM 

TICK-1 #-TXRCK)/TICK) 

TICK«SORTCTICK*TICK) 

IF<TICK-#0005>  570#  570#  400 
570  CONTINUE 

TCKM )*  1 0000# 

P/0  OF  ISTARCX)#  ICR)  AND  TEMPERATURE 
PRINT  3 

PRINTlL#"  ITERATION#  PRINT*1#  NEW  FILE*2 
READiNP 

IFCNP-1)  590#  575#  573 
573  CLOSE  FILE  STOW 
GO  TO  10tl 

575  PRINT*"  R  ISTARCX)  ICR) 

DO  580  K=1 #  KM 

580  PRINT  1#  RCK),  TISCK)#  TICK)#  TCK) 

PRINT  3 

THICK  INVERSION 

ABSORPTION  COEFFICIENT  AND  LOCAL  ISTARCX) 
590  DIT«=0. 

DI  Sa0  • 

DO  600  K»1 #  KM 

600  ABSCK)*TICK)/CCW5/CEXPCC2V/TCK>  >-l . ) > 
TISCKM)*0. 

DO  800  Ksl#  KMX 

KXlaK+1 

XaRCK) 

CALL  INT3PC  TI #  DTI S) 

TISCK)»0. 

DO  640  KX=K#KMX 
640  TISCK)«TISCK)42.*DTISCKX) 


CALL  INT3PC ABS#DTXR) 
TXRCK>=0. 


TEMP" 
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01010 

01020 

01030 

01040 

01050 

01060 

01070 

01080 

01090 

01100 

01110 

01120 

01130 

01140 

01150 

01160 

01170 

01180 

01190 

01200 

01210 

01220 

01230 

01240 

01250 


00010 

00020 

00030 

00040 

00050 


DO  660  KX«K#KMX 
660  TXR<KX*1 >»TXR<KX)+DTXRCKX> 

COMPUTED  0 THICK  I(X)  AND  NEW  ISTARCX) 

TIA<K)e0. 

DO  680  KX«K*KMX 

DTIA*DT1S<KX)*<EXPC-TXR(KM)+TXRCKX) >+EXP(-TXR<KM)-TXRCKX> > > 
680  TIA(K)gsTIA(K)40TIA 

TIS<K>sTIS<K)*TIXCK>/TIA<K> 

*  COMPARISON  OF  MEASURED  AND  COMPUTED  OTHICK  IfX> 

DI*TIA<K>-TIX<K> 

DIT=DIT+DI 
DI S»DI S*DI**2 
800  CONTINUE 

DITaDIT/KMX 
DIS*SORT<DI S/KMX) 

PRINTt"  AVG  DIF  SIGMA- 

PRINT  1*  DIT#  DIS 
CALL  FITXCKZ*N,NNZ> 

CONTINUE 

CALL  IOFX(KZ#N»NNZ) 

GO  TO  220 
1000  STOP 
...  END 


FUNCTION  ROOT<U#V> 
REAL  ROOT 

ROOT= SORT  <U**2-V**2> 

RETURN 

END 
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00010 

00000 

00030 

00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

00160 

00170 

00180 

00190 


00010 
00020 
00030 
00040 
00050 
00060 
00070 
00080 
00090 
00100 
001  10 
00120 
00130 
00140 
00150 
00160 
00170 
00180 
00190 
00200 
00210 
00220 
00230 
00240 
00250 
00260 
00270 
00280 


SUBROUTINE  10FXCKZ*  N#  NNZ> 

COMMON  ARC 6*5) *  MZC5)#  CC6*2)#  RC5D#  TIS<51>#  TIX<51> 

COMMON  ABS(51 )#  T(51)»  TIC51)#  TIAC5D 

COMMON  TXRC51  >  #  DTISC5D#  DTXRC51)#  DR#X#K#KX1 #KMX#  KM 

COMMON  PHIC17*6>#AMC6*6)#C0FC6#7)#YC17)#ABC6#1 > 

KFsMZCI) 

DO  140  K»l#  KF 
TISCK)=ARC 1#1) 

DO  140  J*2#  N 

140  TISCK)*TISCK)+ARC J,1  )*RCK)  **C2*C  J- 1  )  ) 

DO  165  NZ*2*  NNZ 

KI*MZCNZ)-KZ+1 

KF-MZCNZ) 

DO  165  K=KI  #  KF 

150  TISCK)*ARC 1#NZ> 

DO  165  J=2#N 

1 65  TISCK)*TISCK)+ARC J#NZ)*CRCK)**C J-l )) 

RETURN 

END 


SUBROUTINE  THINCN#NNZ) 

COMMON  ARC  6#  5)  #  MZC  5)*  CC6*2>#  RC51 )#  TI SCSI)#  TIXC51) 
COMMON  ABSC51).  TC 51 ) #  TIC51 )#  TIAC51) 

COMMON  TXRC51 ) •  DTISC 51 )#  DTXRC51 >  #  DR»X#K#KX1 #KMX*KM 
COMMON  PHI C 1 7#  6) #  AMC  6#  6) #COFC  6#  7) #  YC 17>#  ABC  6# 1 ) 

OTHIN  ICX)  TO  ICR) 

DO  560  K*  1  #  KMX 
TICK)b0. 

NZ*0 

340  NZ»NZ*1 

X»RCMZCNZ)) 

IFCX-RCK))  340#  340#  350 
350  DEN«RCK) 

LIM  =  3 

IFCNZ-1)  1000#  400#  500 
360  DEN** 1 

370  TICK)*TI CK)+TI Z 
NZ=NZ+1 
LIM=  1 

I FCNNZ“NZ) 560#  500#  500 
380  TICK)=TICK)-TIZ 
X=RCMZCNZ ) ) 

LIM*2 
GO  TO  500 

400  TIZ*2»  *ARC  2# 1 ) 

DO  420  J»3#  N 
*  CP=CCJ,1) 

TIP=CP*C  X**C2+ J- 4) > 
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00290 
00300 
00310 
00320 
00330 
00340 
00350 
00360 
00370 
00380 
00390 
00400 
00410 
00420 
00430 
00  440 
00450 
00460 
00470 
00480 
00490 
00500 
00510 
00520 
00530 
00540 
00550 


00010 

00020 

00030 

00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

00160 


DO  410  JP=3#  J 
CP«CP*CCCJ-JP+2>#  1 ) 

410  TIP=TIP+CP*(X+*C2*( J- JP> )  )♦<  RCK)**C2*( JP-2) > ) 
420  TIZ»TIZ*TIP*ARC J#  1 1 
TIZ»TIZ*ROOTCX#ROO> 

GO  TO  360 

500  TIZsCC3#2>*ARC3#NZ>+CC4#2>*AR<4#NZ>*X 
IFCN-4)  1000#  540#  510 
510  DO  530  J=5#  N 
CP=C< J#2) 

TIPsCP*(X**(J-3>) 

DO  520  JP=5#  J#  2 
CP=CP*CCCJ-JP+3>#2> 

520  TIP»TIP+CP*<X**(J-JP>)*<RCK)**( JP-3>> 

530  TIZeTIZ+TIP*ARCJ#NZ> 

540  TlZ=TIZ*ROOTCX#RCK>> 

TILG»AR( 2#NZ> 

CLG»t 

DO  550  J=4#  N#  2 
CLG=CCJ#2>*CLG 

550  TILG=TILG+CLG*CRCK)**C J-2>  >*ARC J#NZ) 

TIZ=TIZ*TILG*ALOGCCX+ROOTCX#RCK>>>/DEN> 

GO  TO  C 380#  370#  360>#  LIM 
560  CONTINUE 
RETURN 
1000  STOP 
END 


SUBROUTINE  POLYTEM 

COMMON  ARC  6# 5)#  MZC5>#  CC6#2>#  R( 51 >  »  TISC51 )»  TIXC51) 
COMMON  ABSC51)#  TC51 )#  TIC51 >#  TIAC51) 

COMMON  TXRC51 >#  DTISC51 >#  DTXRC51 >»  DR#X»K#KXl »KMX#KM 
COMMON  PHIC17#6>#AMC6#6>#COFC6#7>#YC17)#ABC6#l > 

COMMON  AJC7>#  JORD#JM 
TPW=TCK) 

TXRCK>=A JC  1  > 

DTXRCK)*AJC2) 

DO  100  J*2#  JORD 
DTXR(K)=DTXRCK)+TPW*J*A JCJ+f ) 

TXR(K)eTXRCK)+TPW*AJC J) 

100  TPW*TPW*T<K> 

TXR<K>sTXRCK>*TPW*AJCJM> 

RETURN 

END 
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00010 
00020 
00030 
00040 
00050 
00060 
00070 
00080  & 
00090  6 
00100 
00110 
00120 
00130 
00140 
00 150 
00160 
00170 
00180 
00190 
00200 
00210 
00220 
00230 
00240 
002  50 
00260 
00270 
00280 
00290 
00300 
00310 
00320  6 
00330 
00340 
00350 
00360 
00370 
00380 
00390 
00400 


SUBROUTINE  INT3P<Q*DT> 

COMMON  AR<6&5>*  MZ<5>*  CC6*2>*  RC51>*  TIS<51>*  TIXC 51 > 
COMMON  ABS< 51  )*  T<51>*  TI(51>*  TIA<5t> 

COMMON  TXRC51)*  DTISC51)*  DTXRC5D*  DR*X*K»KX1*KMX*KM 
COMMON  PHI<17*.6>*AM<6*6>*COF(6*7>*Y(17>*AB(6*  1> 

DIMENSION  0(51>*  0IC5D*  E(3>*  ERC3>*  DT(51> 

PTSF <  EU*  EL*  EUR  *  ELR*  FU*  FL  *  PH2*  PH  3 > =FU*EUR-FL*ELR 

-CPH3/2)*CEU*EUR-EL*ELR-AL0G<CEU+EUR>/(EL+ELR>>> 

-CPH2/3>*<EUR**3-ELR**3> 

DO  10  KXsK*KM 
10  OKKX)sO(KX) 

IF(K-l)  15*  15*  20 
15  DTCl>«(DR/3>*<2*QI<l)+9I<2>> 

KX2*2 
GO  TO  60 
20  DE&DR/X 

9I(K-1>«0CK-1> 

E(  1  )  **  1 
ER( 1 )»0* 

DO  30  J=2*  3 
ECJ>=E<  J-D+DE 
30  ER( J)sROOT(E( J)*E( 1 ) ) 

PH1s(0I<K+1)-0I(K-1))/(2*DE) 

PH2*<9I<K-M)-2*0I(K>+9I<K-1>>/DE**2 

PH3=PH1-E(1)*PH2 

DT<K)«X*PTSF(EC2>*E(1 )*ER(2)*ER(1 >*0I(K+1),0I(K>*PH2,PH3) 
IF(KMX-K)  70*  70*  40 
40  PHl»C9ICK+2>-QI<K>)/C2*DE) 

PH2=(QI(K+2)-2»0l(K+l)+0I(K))/DE**2 
PH3*PH 1 -EC  2 ) *PH2 

DT(KX1>*X*PTSFCEC3>*E(2>*ER(3>*ERC2>*0ICK+2>*Q1CK+1>* 

PH2*PH3) 

IFCKMX-KXl)  70*  70*  50 
50  DO  55  KX=KX1*  KM 

55  OI(KX>=OI(KX)*R(KX>/ROOT(R<KX)*X) 

KX2-KX1+1 

60  DO  65  KX*KX2*  KMX 

65  DT(KX>a<DR/12>*(-0I(KX-l >+8*0I <KX)+5*9I CKX+1 >> 

70  RETURN 
END 
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00010 

V  ' 

.  1 

SUBROUTINE  FITXCKZ*N#NNZ > 

00020 

COMMON  ARC 6# 5) *  MZC5>  #  CC6#2>#  RC51 >#  TI SC  51 >#  TIXC51 > 

00030 

COMMON  ABSC5D#  TC51>#  TIC51>#  TIAC51 > 

00040 

COMMON  TXRC51 >#  DTISC51)#  DTXRC51 >  *  DR#X#K#KX1#KMX#KM 

00050 

COMMON  PHI C 1 7#  6) # AMC  6#  6) #C0FC6#  7># YC 17>#ABC6#1> 

00060 

NB*1 

00070 

NA»N 

00080 

DO  950  NZ*1#  NNZ 

00090 

IFCNZ-1 >  1000*  400 «  500 

00100 

400 

MS*0 

001 10 

M=0 

00120 

GO  TO  650 

00130 

500 

MSsMZCNZ-1 )-4 

00140 

IFCMS)  550#  600#  600 

00150 

550 

MS=0 

00160 

M*MZCNZ-1)-1 

00170 

GO  TO  650 

00180 

600 

M*3 

00190 

650 

IFCKM-CMZCNZ>+3>>  700#  750#  750 

00200 

700 

M=M+CKZ+1>*CKM-MZCNZ>> 

00210 

GO  TO  800 

00220 

750 

MsM+CKZ+1 )+3 

00230  . 

.800 

REF»DR*CM-1 ) 

00240 

DO  880  1*1  #M 

00250 

PHI C I # 1 )* 1  • 

00260 

YC I )=TISC I  ♦MS)  , 

00270 

IFCNZ-1 >  1000#  820#  840 

00280  . 

820 

DO  830  J*2#N 

00290 

830 

PHICI. J)sCRCI)/REF)**C2*CJ-l >> 

00300 

GO  TO  880 

00310 

840 

DO  850  J=2#N 

00320 

850 

PHICI#J)sCDR*CI-l>/REF)**CJ-l> 

00330 

880 

continue 

00340 

890 

DO  892  JL*1#  N 

00350 

DO  892  JK=1#N 

00360 

AMC JK# JL)*0» 

00370 

DO  892  1*1 #M 

00380 

892 

AMCJK#  JL)*AMC  JK#  JL1+PHICI#  JL>*PHICI#  JK) 

00390 

DO  896  JK* 1 #N 

00400 

ABC JK#1 )*0. 

00410 

DO  896  I* 1 #M 

00420 

896 

ABCJK#1 )*ABCJK#1 )»YC I >*PHI C I, JK) 

00430 

CALL  LINE0CAM#AB#NA#NB#6> 

00440 

DO  915  J*  2#N 

00450 

IFCNZ-1)  1000#  900#  910 

00460 

900 

ABC J# 1 )*ABC J# 1 )*C  REF**C2*C 1- J) ) ) 

00470 

GO  TO  915 

00480 

910 

ABC J# 1 )*ABC J# 1 ) *  C  REF**  C 1 - J> ) 

97 


00490 

00500 

00510 

00520 

00530 

00540 

00550 

00560 

00570 

00580 

00590 

00600 

00610 

00620 

00630 

00640 

00650 

00660 

00670 


915  CONTINUE 

IF<MS>  1000«  9424  920 
920  DO  930  Jsl«N 
JP=J-1 

CALL  BICOFC04 JP4COF) 

DO  930  JR=l4  J 
930  COF<  JR4 J+l >=COFC  JR4 1 > 

DO  940  JRsl,  N 
ARCJR4NZ>=0. 

DO  940  JsJR*  N 

940  ARCJR4NZ>=AR(JR4NZ>+((-| >**( J- JR))*AB( J4 1 )*(RtM 
6  ZCNZ-l )-3)**CJ-JR>>*C0FC JR4( J+! )) 

GO  TO  950 
942  DO  946  J**l4N 
946  ARCJ4NZ>sABCJ41) 

9S0  CONTINUE 
RETURN 
1000  STOP 
END 
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